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Chapter  1 

INTRODUCTION 


1.1  Problem  Description 

The  purpose  of  this  report  is  to  develop  using  SIMULINK  an  aircraft  simulation  that 
couples  the  nonlinear  equations  of  motion  with  the  nonlinear  aerodynamic  and  engine 
performance  data.  The  simulation  has  been  designed  using  data  from  actual  flight 
testing  of  the  McDonnell  Douglas  F-15  Eagle.  The  ultimate  objective  is  to  provide  a 
generic  framework  for  the  development  of  nonlinear  simulations  of  other  aircraft  with 
minimal  required  changes.  Design  of  a  longitudinal  autopilot  using  the  Total  Energy 
Control  System  (TECS)  concept  is  accomplished  using  the  linearized  model  and  later 
validated  with  the  nonlinear  model. 

1.2  Analysis  of  Data  Provided 

The  following  sections  summarize  the  model  characteristics  and  the  relevant  support 
modules  for  the  F-15  described  in  [5].  The  numerical  data  were  provided  in  the  format 
of  a  Genesis  simulation  used  at  Wright-Patterson  and  the  coding  was  subsequently 
converted  from  FORTRAN  to  MATLAB.  A  detailed  description  of  the  model  follows 
in  Chapter  2. 

The  evaluation  of  the  open-loop  nonlinear  model  and  the  design  of  the  TECS 
control  law  will  be  accomplished  at  the  two  flight  points  shown  in  Table  1.1.  Other 
flight  conditions  throughout  the  flight  envelope  will  also  be  investigated. 

1.2.1  Model  Characteristics 

The  model  is  an  integration  of  several  modules,  each  performing  a  specific  function. 
These  modules  include  the  aerodynamic,  propulsion,  and  atmospheric  models  (see 
Appendix  E),  as  well  as  the  nonlinear  equations  of  motion.  Modeling  of  the  control 
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Table  1.1;  Selected  Flight  Conditions 


Flight 

Altitude 

Vtas 

Mach 

Point 

(ft) 

(ft/s) 

Number 

1 

9,800 

539.1 

0.5 

2 

30,000 

497.3 

0.5 

surface  and  actuator  dynamics  as  well  as  any  sensor  dynamics  is  performed  where 
necessary  in  the  SIMULINK  environment.  The  closed-loop  control  law,  in  this  case 
TECS,  is  also  implemented  in  SIMULINK.  Figure  1.1  summarizes  the  integration  of 
these  components  to  form  the  complete  system  model. 


Figure  1.1:  Integration  of  System  Model  Components. 


The  aircraft  modeled  is  the  McDonnell  Douglas  F-15  Eagle,  the  state-of-the-art 
in  current-day  operational,  high-performance  fighter  aircraft.  It  is  powered  by  two 
afterburning  turbofan  engines,  each  providing  approximately  32,000  lb  of  thrust  at 
maximum  power.  The  primary  flight  control  surfaces  include  horizontal  stabilators 
capable  of  both  symmetric  and  differential  movement,  conventional  ailerons,  and  twin 
vertical  rudders.  There  are  a  total  of  six  actuators-two  stabilators,  two  ailerons,  and 
two  rudders.  All  actuators  are  modeled  identically  with  rate  limits  of  24  deg/sec  and 
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first-order  response  characteristics  of 


G(.)- 


20 

s  -\-  20 


The  individual  surface  position  limits  and  sign  conventions  for  positive  deflection 
are  summarized  in  Table  1.2.  The  aircraft  mass  and  geometry  characteristics  are 
summarized  in  Table  1.3. 


Table  1.2:  Flight  Control  Sur{a,ces  Characteristics 


Control  Surface 

Symbol 

Limits 

Sign  Convention  (±) 

Symmetric  stabilator 

Sh 

±20° 

Trailing  edge  down 

Differential  stabilator 

±15°/  -  25° 

Left  trailing  edge  down 

Aileron 

Sa 

±20° 

Left  trailing  edge  down 

Rudder 

±30° 

Trailing  edge  left 

Table  1.3:  F-15  Mass  and  Geometry  Characteristics 


Parameter 

Symbol 

Units 

Value 

Wing  area 

s 

ft^ 

608.0 

Wing  span 

b 

ft 

42.8 

Mean  aerodynamic  chord 

c 

ft 

15.95 

Aircraft  weight 

W 

lb 

45,000.0 

Moments  of  inertia 

h 

slug  —  ft^ 

28,700.0 

ly 

slug  —  ft^ 

165,100.0 

h 

slug  —  ft^ 

187,900.0 

Products  of  inertia 

^XZ 

slug  —  ft^ 

-520.0 

Jxy 

slug  —  ft^ 

0.0 

Jyz 

slug  —  ft^ 

0.0 
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1.2.2  Aerodynamic  Model 

The  aircraft  aerodynamics  are  modeled  using  a  combination  of  multidimensional  ta¬ 
bles  and  linear  interpolation  to  form  nonlinear  function  generators.  The  highly  non¬ 
linear  aerodynamics  encountered  in  the  extreme  portions  of  the  aircraft  flight  envelope 
are  therefore  represented,  which  will  be  valuable  in  validating  the  control  law  after  it 
has  met  the  requirements  for  the  linearized  model.  Most  of  the  aerodynamic  quanti¬ 
ties  are  a  function  of  Mach  number  M,  and  some  combination  of  angle  of  attack  a, 
sideslip  angle  /?,  and  symmetric  stabilator  deflection  Sh- 

The  aerodynamic  model  calculates  the  nondimensional  force  and  moment  coeffi¬ 
cients,  which  are  then  used  to  calculate  the  total  associated  forces  and  moments.  The 
equations  used  for  the  coefficients  are 


•  Coefficients  of  forces 

Cl  =  Clo  +  ^Chn^nz 

Cd  =  Cdo  +  l^COal,  +  ^Conoz 

Cy  =  Cyo  +  Cys^ Sa  +  Cy^^ do  -  ACyc^Ksh^ 


•  Coefficients  of  moments 

Ce  —  Ce^  +  Cif^SA  +  Ceg^So- ACee^Ks^^  + ^{CepP  +  Ctrr) 

Cm  —  Cmo  A  ACmn^'lT'z  ’^{Cmgfl  A  Cma^  A  Cl^ANo) 

Cn  —  Cno  +  Cnsj^^A  A  C„f^  Sj)  -f  ACn^^KsR^  +  2v(^«pP  + 

The  terms  C,  AC,  AN,  and  K  are  outputs  from  the  function  generation  routines, 
and  are  either  calculated  directly  or  by  linear  interpolation  of  the  tabular  data.  In 
the  case  of  the  F-15,  all  coefficients  are  determined  from  tables  except  AC'£)„„^  and 
Cd  for  Oi  >  40°,  which  are  calculated  directly.  The  parameters  affecting  each  of  the 
coefficients  are  summarized  in  Table  1.4.  The  total  forces  and  moments  are  calculated 
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Table  1.4:  Independent  Parameters  Affecting  the  Aerodynamic  Coefficients 
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from  the  equations 

L  =  qSCL 
D  =  qSCo 
Y  -  qSCv 

SL  -  qSbCi 
EM  =  qSWra 
SiV  -  qSbCn 

where  q  —  is  the  dynamic  pressure,  S  is  the  wing  area,  b  is  the  wing  span,  and 
c  is  the  mean  aerodynamic  chord. 

1.5.3  Propulsion  Model 

The  propulsion  system  model  consists  of  two  distinct  engine  models.  Although  the 
two  engines  are  similar,  they  are  not  exactly  identical — i.e.  for  a  given  throttle  setting, 
the  thrust  produced  may  vary  slightly.  The  engine  thrust  vectors  are  aligned  with 
the  aircraft  x  body-axis,  and  the  thrust  produced  is  a  function  of  altitude  h,  Mach 
number  M,  and  throttle  setting  Spla-  Each  engine  is  modeled  as  a  nonlinear  system 
having  two  separate  sections — a  core  engine  and  an  afterburner  (augmentor),  each 
with  its  associated  sequencing  logic. 

The  throttle  position  inputs  to  the  engine  model  are  in  degrees  of  power- level- angle 
(PLA),  with  a  minimum  position  of  20®  and  a  maximum  of  127°.  The  core  section 
responds  to  settings  up  to  83°,  while  the  afterburner  section  begins  to  respond  at  a 
position  of  91°.  The  core  model  has  first-order  dynamics  and  rate  limiting  to  model 
spool-up  effects,  while  the  afterburner  has  a  rate  limiter  and  sequencing  logic  to  model 
the  fuel  pump  and  pressure  regulator  effects. 

1.2.4  Atmospheric  Model 

The  atmospheric  model’s  data  is  based  on  tables  from  the  U.S.  Standard  Atmosphere 
(1962).  This  model  calculates  values  for  the  speed  of  sound,  the  acceleration  due  to 
gravity,  air  density,  viscosity,  and  ambient  static  pressure  and  temperature  based  on 
the  aircraft  altitude.  Linear  interpolation  is  used  between  table  values  for  altitudes 
from  0  to  90  km. 
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1.2.5  Equations  of  Motion 

The  nonlinear  equations  of  motion  used  in  the  system  model  are  based  on  the  deriva¬ 
tion  by  Duke,  Antoniewicz,  and  Krambeer  in  [6j.  These  equations  model  the  six- 
degree-of-freedom  dynamics  of  a  rigid  aircraft  flying  over  a  flat,  nonrotating  Earth. 
The  derivations  are  outlined  in  Chapter  2. 


Chapter  2 

THE  NONLINEAR  F-15  MODEL 


Linear  aircraft  models  are  useful  in  the  early  development  of  control-law  design 
and  in  the  analysis  of  vehicle  dynamics  over  many  flight  conditions.  However,  because 
such  models  are  only  approximations  of  the  aircraft  behavior,  it  is  valuable  to  the 
engineer  to  be  able  to  verify  the  performance  of  a  control  law  with  the  corresponding 
nonlinear  model. 

2.1  Derivation  of  Nonlinear  State  Equations 

Most  linearized  aircraft  models  do  not  necessarily  decouple  into  the  longitudinal  and 
lateral  modes.  However,  if  we  include  additional  simplifying  assumptions  such  as 
vehicle  symmetry  or  a  specific  reference  trajectory  such  as  straight-and-level  flight, 
these  modes  can  be  decoupled.  Motion  of  an  aircraft  can  be  modeled  by  the  responses 
of  the  rigid-body  dynamics  described  by  a  set  of  six  nonlinear  simultaneous  second- 
order  differential  equations.  In  the  model  to  be  developed  here,  these  equations  will 
be  used  without  assuming  any  reference  trajectory  or  vehicle  symmetry.  Instead,  the 
model  will  assume  a  rigid  aircraft  of  constant  mass  flying  over  a  flat,  nonrotating 
earth. 

8.1.1  Reference  Systems 

The  primary  reference  systems  are  the  body-,  the  wind-,  and  the  vehicle- carried, 
vertical-axis  systems.  Each  has  its  advantages  for  use  with  a  particular  set  of  equa¬ 
tions  or  variables. 

The  rotational  equations  of  motion  are  most  easily  referenced  to  the  body-axis 
system.  The  body-axis  rotational  rates  are  measureable  within  the  aircraft  by  sensors 
fixed  in  the  body  frame.  The  body-axis  system  has  its  origin  at  the  aircraft  center  of 
gravity,  with  the  a;-axis  directed  out  the  aircraft  nose,  the  ^^axis  out  the  right  wing, 
and  the  ^r-axis  out  the  bottom  of  the  aircraft.  The  positive  direction  for  the  body-axis 
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rates  [p,  q  and  r),  velocities  (w,  v  and  w),  and  moments  (L,  Mand  Pf)  are  shown  in 
Figure  2.1. 


The  wind-axis  system  is  used  for  the  translational  equations  of  motion.  Because 
aerodynamic  forces  act  in  the  direction  of  the  wind  axes,  such  quantities  as  the  angle 
of  attack  a,  the  total  velocity  V,  and  the  sideslip  angle  /S  are  directly  measurable, 
or  closely  related  to  directly  measureable  quantities,  in  the  aircraft.  The  x-axis  in 
the  wind-axis  system  is  aligned  with  the  aircraft  velocity  vector,  with  the  y  and  z 
axes  exiting  the  right  side  and  the  bottom  of  the  aircraft  respectively.  Because  both 
the  wind-  and  body-axis  systems  have  their  origin  at  the  center  of  gravity,  their 
orientation  can  be  defined  from  ^  and  a.  Components  of  the  total  velocity  V  can  be 
expressed  in  terms  of  the  body-axis  velocities  as, 

V  cos  a  cos  (0 

Vsm^  (2-1) 

V  sin  a  cos 

The  vehicle- carried,  vertical-axis  system  also  has  its  origin  at  the  aircraft  center  of 
gravity.  This  reference  system  is  primarily  useful  for  its  orientation  to  the  body-axis 
system,  which  is  defined  by  the  Euler  angles  ?/>,  9  and  ^  (Figure  2.2).  The  vehicle- 
carried,  vertical-axis  system  is  defined  by  having  the  a:-axis  directed  north,  the  j^axis 
directed  east,  and  the  z-axis  directed  down.  This  axis  system  is  simply  the  earth-fixed 
reference  system  translated  to  the  aircraft  center  of  gravity. 
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*2’  *6 


Figure  2.2:  Orientation  of  Vehicle- Carried,  Vertical-Axis  System  to  the  Body-Axis 
System. 

Based  on  these  axis  systems,  the  aircraft  dynamics  can  be  described  by  12  states, 
which  will  be  divided  into  four  sets  of  three  variables,  representing  the  vehicle  ro¬ 
tational  velocity  (p,  q,  and  r),  the  vehicle  translational  velocity  (V,  a,  and  /?),  the 
vehicle  attitude  ((/!>,  6,  and  ^),  and  the  vehicle  position  [h,  x,  and  y). 

2.1.2  Rotational  Accelerations 

The  rotational  acceleration  terms  p,  q  and  f  are  derived  from  the  moment  equation 

M  =  (2.2) 

where  M  is  the  total  moment  on  the  aircraft  and  H  is  the  total  angular  momentum 
of  the  aircraft.  For  a  moving  reference  frame  such  as  the  body-axis  system,  the  time 
derivative  operator  6 /St  is  used  in  place  of  d/dt,  and  the  total  angular  momentum  is 
given  by  H  =  /fi,  where  I  is  the  inertia  tensor  and  Q  is  the  rotational  velocity  vector. 
The  result  is 


M  =  +  nx  (in) 

St 


(2.3) 


where 
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'  ST 

L  -j-  Lj' 

M  = 

SM 

M  -f-  Mj' 

SiV 

N  +  Nt 

(2.4) 


with  L,  Mand  N  defined  as  the  aerodynamic  total  moments  about  the  x,  y  and  z  body 
axes  respectively,  and  Lt-,  Mt  and  Nt  defined  as  the  sums  of  all  the  thrust-induced 
moments.  The  inertia  tensor  /  is  defined  as 


7  = 


-I. 


xy 


txy 


lyz 


‘yz 


(2.5) 


where  the  moments  of  inertia  about  the  x,  y  and  2;  body  axes  and  the  products  of 
inertia  in  the  x-y,  x-z  and  y-z  body-axis  planes  are 


•  Moments  of  inertia: 


Products  of  inertia: 


4  =  /  {y^  +  dm 

Jb 

ly  =  I  +  z^)  dm 

J  B 

4  =  /  +  y^)dm 

J  B 

Ixy  —  I  xy  dm 

JB 

Ixz  —  I  xz  dm 

Jb 

lyz  =  /  yzdm 

Jb 


The  term  is  a  vector  of  the  rotational  rates  p,  q  and  r  about  the  x,  y  and  2 
body  axes  respectively.  Because  constant  mass  is  assumed  for  the  aircraft,  the  inertia 
tensor  is  constant  with  respect  to  time,  and  equation  (2.3)  can  be  solved  for  the 
rotational  accelerations, 


12 


[M  -  X  (/fi)] 


(2.6) 


In  order  to  simplify  the  expansion  of  this  expression,  the  inverse  of  the  inertia  tensor 
will  be  defined  as 


where 


/-i 


1 

det  I 


h  h  h 
h  h  h 
h  h  h 


det  I  —  Ixlylz  “^lyzlxzlxy 

11  =  lylz  lyz 

12  —  Ixy^z  "b  iyzlxz 
Iz  —  Ixylyz  “b  lyJxz 
h  -  Uz  -  4 

-^5  —  ^x^yz  ^xy^xz 


•^6  —  ^x  ^y  ^ 1 


2 

xy 


(2.7) 


and  the  terms  Dy  and  Dz  are  defined  as 

Dx  —  Iz  ly 

Dy  —  Ix  Iz 
D z  —  ^y  Ix 

Expressions  for  the  rotational  accelerations  in  equation  (2.6)  can  now  be  expanded 
into  the  following  set  of  scalar  equations, 
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'  p  =  g^[(SL)/i  +  (SM)/2  +  (SiV)/3-p^(4./2- Wa) 

~\~PQ(^Ixz^1  Iyz^2  H” 

i,^yz^l  Ixy^s)  9^(-Dj;/i  Ixyl2  4"  ^ar^^s)  ^  ^xzl2)\ 

q  -  ^[(EZ)/2  +  (SM)/4  +  (SAr)/5-p2(4^^^ 

^  ~\~PQ{lxzh  “  -^yJK^4  ”  £^.2^5)  ~  P'i^il^xyh  +  Dyl^  —  lyzh)  (2*8) 

4~9  (^y2^2  /ry^s)  q'f'(^Dxl2  ^ci/^4  H“  Ixz^S^  ^  (^3/2^2  -^r2-^4)] 

r  -  5^[(SL)/3  +  (EM)/5  +  (SiV)/6-p2(4^^^ 

~\~pq(^Ixz^3  ^yz^S  -^2-4)  P^{4y-4  H“  ^y^S  ^yz^6^ 

.  4"9  -4y4)  -4y4  "f"  424)  ^  (4^4  424)] 

2.L3  Translational  Accelerations 

Derivation  of  the  translational  acceleration  terms  is  based  on  the  force  equation 

F  =  I  (mV)  (2.9) 

where  F  represents  the  total  force  acting  on  the  aircraft,  m  is  the  aircraft  mass,  and 
V  is  the  total  velocity.  By  assuming  m  to  be  constant  and  using  the  time  derivative 
operator  for  a  rotating  reference  frame  6 /St,  this  expression  can  be  expanded  to 

C 

F=xm(— V  +  OxV)  (2.10) 

where  F  =  [SX  SF  SZ]^,  a  vector  containing  the  sums  of  the  aerodynamic,  thrust, 
and  gravitational  forces  in  the  x,  y,  and  .^body  axes  respectively,  and  V  ==  [u  u  u;]^. 
With  some  manipulation,  equation  (2.10)  can  be  solved  for  the  body-axis  translational 
accelerations. 


— F-Q  X  V 
m 


(2.11) 


However,  as  previously  mentioned,  the  desired  form  for  the  translational  acceler¬ 
ations  is  in  the  wind-axis  system — in  terms  of  the  total  velocity  magnitude  V,  the 
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angle  of  attack  a,  and  the  sideslip  angle  p.  The  relationships  between  the  trans¬ 
lational  accelerations  of  the  body-  and  wind-axis  systems  are  shown  in  equations 
(2.1).  Derivation  of  the  scalar  equations  for  the  wind-axis  translational  accelerations 
is  involved  and  will  not  be  shown  here.  The  resulting  equations  are 

V  —  ^[— Dcos /3  +  F  sin -f  cosa  cos /? -f- Irsin  jd  +  Zy  sina  cos /? 

—mg{cos  a  cos  /9  sin  0  —  sin  /?  sin  (j>  cos  ^  —  sin  a  cos  /S  cos  ^  cos  0)] 

a  =  Z,  +  Zt  cosof  —  Xy  sina -f- m^(cosQ;cos  (^cos  0  +  sin  q;  sin0)] 

<  —  tan /9(p  cos  a -f  r  sin  a)  (2-12) 

^  =  -^[D  sin  /9  -f-  F  cos  ^  —  Xt  cos  a  sin  ^  -f-  It  cos  /3  —  Zt  sin  a  sin  /3 

-|-m5r(cos  a  sin  sin  0  -|-  cos  ^  sin  <j)  cos  ^  —  sin  a  sin  ^  cos  (j)  cos  0)] 
d-p  sin  a  —  r  cos  a 

where  D  is  the  total  aerodynamic  drag,  F  is  the  total  aerodynamic  side  force,  L  is 
the  total  aerodynamic  lift,  and  Xt,  Yt  and  Zt  are  the  total  thrust  forces  in  the  x,  y 
and  2r  body-axis  directions  respectively. 

2.1.4  Attitude  Rates 

The  Euler  angle  rates  in  the  earth-fixed  axis  system  and  the  rotational  velocities  in 
the  body-axis  system  are  related  by  a  transformation  matrix  T,  where 

1  0  —  sin  ^ 

T  =  0  cos4>  sin  ^  cos  ^  (2-13) 

0  —  sin  (j)  cos  (f)  cos  9 

The  equation  for  the  transformation  from  the  earth-fixed  system  to  the  body-axis 
system  is 

fi  =  r(4E)  (2.14) 

with  E  representing  a  vector  of  the  Euler  angles,  E  =  6  By  a  simple 

rearrangement,  equations  for  the  attitude  rates  are 


15 


lE  =  T-‘n  (2.15) 

which  can  be  expanded  into  the  following  scalar  equations 

^  —  p qsin<f)t&n6  +  r  cos  (f>ta,n9 
6  —  q  cos  <f>  —  r  sin  ^ 

Ip  —  q  sin  <f>  sec  9  r  cos  (p  sec  9 

2.1.5  Earth-Relative  Velocity 
The  earth-relative  velocities  are  related  to  the  body-axis  velocities  by 

V  =  Lbv(^K)  (2.17) 

where  R  is  the  location  of  the  aircraft  in  the  earth-axis  system,  R  =  [a;  y  z]^,  and 
z  =  —h.  Note  that  V  is  the  total  velocity  vector  with  components  u,  v  and  w  in  the 
body-axis  system.  The  transformation  from  the  earth-axis  system  to  the  body-axis 
system  is  accomplished  by  the  matrix  Lbv-,  where 


(2.18) 

Equation  (2.17)  is  easily  solved  for  the  earth- relative  velocities  in  terms  of  the  body- 
axis  velocities: 

jK=Ll\Y  (2.19) 

Using  this  equation  and  equation  (2.1),  the  earth-relative  velocities  can  be  solved  in 
terms  of  V,  a  and  (i.  The  resulting  equations  are 


cos  ip  —  sin  Ip  0 
sin  Ip  cos  xp  0 

0  0  1 


cos  9  0  sin  0  10  0 

0  10  0  cos  <p  —  sin 

—  sin  0  0  cos  0  0  sin  cos  <p 


(2.16) 


16 


V (cos  a  cos  /?  sin  0  —  sin  sin  ^  cos  0  —  sin  o;  cos  (3  cos  (f)  cos  9) 

F[cos  a  cos  ^  cos  6  cos  ^  +  sin  ^(sin  ^  sin  0  cos  ^  —  cos  <;;!>  sin  0) 

+  sin  a  cos  /?(cos  ^sin  0  cos  »/>  +  sin  sin  ^)]  (2.20) 

^^[cos  a  cos  ^  cos  0  sin  ^  +  sin  /?(cos  <f)  cos  ^  +  sin  <;/>  sin  0  sin  ij;) 

+  sin  a  cos  ,9(cos  (j)  sin  0  sin  i/>  —  sin  ^  cos  ^)] 

2.2  Nonlinear  Simulation  Model 

The  twelve  nonlinear  simultaneous  second-order  diflferential  equations  derived  in  Sec¬ 
tion  2.1  are  implemented  in  the  simulation  in  the  form  of  a  MATLAB  special  function, 
or  s-function.  The  intricacies  of  this  function  are  left  to  the  reader  to  learn;  however, 
the  basic  format  and  flow  will  be  described  here.  Integration  with  the  simulation 
program  SIMULINK  will  also  be  discussed  in  general  terms  to  give  the  reader  an 
overall  picture  of  how  the  various  components  of  the  simulation  interact. 

2.2.1  Component  Integration — the  S-Function 

Systems  saved  under  SIMULINK  as  s-functions  behave  like  MATLAB  m-functions 
and  can  be  called  from  user-written  routines  or  from  the  command  line.  The  s- 
function  makes  it  possible  for  the  user  to  write  customized  routines  for  simulation, 
linearization,  and  parameter  estimation.  As  described  in  Chapter  1,  the  F-15  dy¬ 
namics  are  contained  in  two  primary  routines,  f  25aero  and  f  25eng,  quantifying  the 
aircraft  aerodynamic  and  propulsion  characteristics  respectively.  The  routine  atmos 
calculates  the  current  atmospheric  parameters  such  as  density  and  temperature  to 
be  used  in  these  routines.  Using  the  updated  values  of  the  system  states  and  inputs, 
the  necessary  aerodynamic  and  propulsion  parameters  are  calculated,  which  are  then 
used  in  the  calculation  of  the  twelve  state  derivatives.  All  three  of  these  routines  were 
provided  in  FORTRAN  format  from  the  AIAA  Design  Challenge  [5].  The  s-function 
used  for  the  open-loop  simulation  is  shown  in  Appendix  A. 

Initial  conditions  for  the  aircraft  states  are  also  specified  in  the  s-function.  These 
values  are  of  particular  importance  when  determining  the  aircraft  trim  point  for  a 
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given  flight  condition.  This  will  be  discussed  in  further  detail  in  Section  2.3.2.  The 
s-function  also  allows  the  system  output  values  to  be  specified.  For  an  open-loop 
model,  the  outputs  selected  are  arbitrary.  In  this  case,  the  twelve  states  are  selected 
in  order  to  validate  the  equilibrium  state  responses  at  the  two  flight  conditions. 

2.2.2  The  Simulink  Model 

While  the  s-function  allows  us  to  specify  a  particular  set  of  ordinary  differential 
equations,  the  SIMULINK  environment  allows  the  user  to  design  the  overall  system 
including  other  more  generic  operations  such  as  inputs,  filters,  and  feedback  loops. 

The  open-loop  SIMULINK  system  is  shown  in  Figure  2.3.  The  primary  use  for 
this  simulation  format  is  to  establish  a  trim  point  at  selected  flight  points  and  evaluate 
the  open-loop  aircraft  responses  to  initial  conditions,  control  inputs,  or  disturbances. 
The  initial  conditions  are  specified  in  the  appropriate  s-function,  here  called  f 25sfn 
for  the  open-loop  model.  Additional  inputs  can  be  set  for  the  aircraft.  These  include 
the  aircraft  flight-control  surfaces,  which  consist  of  horizontal  stabilators  capable  of 
symmetric  (DH)  and  differential  (DD)  movement,  conventional  ailerons  (J5A),  and 
twin  vertical  rudders  {DR).  In  this  report,  these  terms  will  be  referred  to  as  8^,  Sd, 
8 A,  and  respectively.  Thrust  from  the  two  engines  is  specified  by  the  right  and 
left  power-level-angles,  PLAR  and  PLAL  respectively,  or  can  be  set  identically  as 
PLASYM.  Again,  the  nomenclature  in  this  text  for  these  variables  is  8plar,  8plal 
and  8pLA-  The  outputs  include  the  twelve  states,  plus  the  angle-of-attack  rate  a. 

2.3  Linearization  of  the  Model 

Design  of  a  feedback  control  system  for  a  nonlinear  system  is  greatly  simplified  by 
the  use  of  a  linearized  model  about  some  nominal  trajectory.  For  aircraft,  it  is 
common  practice  to  use  several  models  linearized  about  various  points  within  the 
flight  envelope  to  design  the  control  system.  Normally,  the  controller  structure  is 
designed  to  be  constant  throughout  the  flight  envelope,  with  some  method  of  gain 
scheduling  used  to  achieve  adequate  performance  at  all  points.  The  linearized  model 
is  therefore  a  significant  tool  in  simplifying  the  development  of  a  satisfactory  controller 
for  the  real-world,  nonlinear  aircraft. 
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Figure  2.3:  SIMULINK  Model  for  Open-Loop  F-15  Simulation. 


2.3.1  Mathematical  Approach 

In  this  section,  the  linearization  of  the  aircraft  equations  of  motion  will  be  summa¬ 
rized.  The  reader  is  referred  to  [6]  for  a  more  detailed  derivation.  Beginning  with 
the  translational  equations  of  motion  (2.12),  by  allowing  the  thrust  and  aerodynamic 
forces  to  be  represented  as  the  total  external  forces  A",  Y  and  Z,  the  result  is 


V  =  — 5'(sin  9  cos  a  cos  /?  —  cos  9  sin  d*  sin  /?  —  cos  9  cos  ^  sin  a  cos 
-j-  —  cos  a  cos  /3  -1-  —  sin  -f-  —  sin  a  cos  B 

d  =  g-pcosatan/?-rsinatan/?-l-5£2s»c«si^±sini^ 

_ X  sing  i  Z  cos  or 

m  V  cos  (3  mV  cos  (3 

$  =  psinof  —  rcosa -f  ■^(sin0cosasin(0-t-cos0sin(^cos^ 

—  cos  9  cos  (j)  sin  a  sin  /?)  —  ^  cos  a  sin  -F  ^  cos  /?  —  ^  sin  a  sin  ^ 

The  linearized  equations  are  now  determined  in  terms  of  the  perturbation  variables 
AF,  Aa,  A/?,  Ap,  A^,  Ar,  AX,  AF,  AZ,  A9  and  A^  in  a  symmetric  climb  with 


(2.21) 


F  =  F,  -F  AF 


(2.22) 
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where  Vo  is  the  aircraft  trim  velocity  (a  constant).  Similar  relationships  as  in  equation 
(2.22)  exist  for  all  the  other  variables  listed  above.  In  this  example,  for  a  symmetric 
climb,  several  trim  values  are  equal  to  zero;  namely,  —  Po  =  (lo  —  '>'0  —  Yg  =  (l)o  — 

=  0.  These  trim  values  may  not  be  zero  for  other  flight  conditions,  such  as  a  level 
turn.  Because  Ah',  Aa,  A/?,  Ap,  Aq,  Ar,  AA,  AY,  AZ,  A9  and  A(f)  are  pertur¬ 
bations  about  their  respective  trim  values,  they  are  always  treated  as  small  quantities. 

By  substituting  the  expanded  variables  of  the  form  of  equation  (2.22)  into  equa¬ 
tions  (2.21),  linearized  equations  of  motion  for  the  perturbed  variables  AF,  Act  and 
A/S  are  derived.  This  is  accomplished  by  assuming  all  second-  and  higher-order 
terms  to  be  negligible  (i.e.  AV  Aa  0)  and  using  the  small  angle  approximations, 
cos  Act  1  and  sin  Act  «  Act.  Leaving  some  relatively  complicated  manipulations  to 
the  reader,  the  equations  become 


AF  =  -pcos(0,-ct,)A0+^AX  +  ^AZ 
<  Act  =  Ag-^sin(0,-ct,)A0-^AA  +  ^AZ  (2.23) 

A/?  =  sin  cto  Ap  —  cos  cto  Ar  +  ^  cos  00  -f  AF 

Further,  by  expressing  the  external  forces  X  and  Z  again  in  terms  of  the  aircraft  lift 
L,  drag  D,  and  thrust  Tand  using  the  perturbation  variables  as  described  in  equation 
(2.22),  the  linearized  equations  of  translational  accelerations  are  reduced  to 


'  AV  =  -g  cos  (0.  -  A9  -  ^  Ao  -  ^  AC  +  Ar 

■  AA  =  (2.24) 


A$  =  sin  cto  Ap  —  cos  Oo  Ar -f  ^  cos  0o  A^  +  Ay 


The  linearized  equations  for  rotational  accelerations  are  found  in  a  similar  manner, 
defining  the  rotational  velocities  (p,  q  and  r)  and  the  moments  (ST,  SM  and  SA^)  in 
terms  of  their  trim  and  perturbation  values.  Using  the  assumed  level  climb  condition, 
we  have  Po  =  =  Tq  =  0  and  STo  —  SMo  =  SA^o  ==  0.  Substituting  the  trim  and 

perturbed  variables  into  equations  (2.8)  and  retaining  only  the  first-order  terms  in 
Ap,  A^  and  Ar,  the  following  linearized  equations  result. 
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A?  -  4j,  Ap  -  4^  Ar  =  A(SM)  (2.25) 

4,Ap-7,,Ar  =  A(SL)  (2.26) 

4.Ar-4,Ap=:  A(SAr)  (2.27) 

By  assuming  a  vehicle  geometry  that  is  symmetric  about  the  body  x-z  plane,  we  have 
4y  —  lyz  —  0  and  equation  (2.25)  can  be  further  reduced  to 

lyy  Aq  =  A(SM)  (2.28) 

Finally,  the  linearized  attitude  (Euler)  rate  equations  are  determined  by  substi¬ 
tuting  into  equations  (2.16)  the  appropriate  trim  and  perturbation  variables.  After 
applying  small  angle  approximations  and  eliminating  higher-order  terms,  the  result¬ 
ing  linearized  equations  are 

'  Ae  =  Aq 

■  =  ^Ar  (2.29) 

A^  =  Ap  +  tan  00  Ar 

2.3.2  Determination  of  the  Trim  Point 

When  a  system  is  nonlinear,  an  operating  point  must  be  chosen  at  which  to  extract 
the  linearized  model.  For  the  F-15  or  any  other  aircraft,  the  equilibrium  (or  trim) 
values  for  the  system  states  and  inputs  are  determined  for  a  given  flight  condition. 
In  this  case,  the  MATLAB  algorithm  trim  is  used.  This  routine  proves  adequate  in 
providing  the  required  data,  however,  several  drawbacks  exist  which  suggest  the  need 
for  an  improved  program  in  future  research. 

The  purpose  of  the  trim  routine  is  to  determine  steady-state  parameters  that 
satisfy  input,  output,  and  state  conditions.  This  is  accomplished  by  searching  for  the 
inputs  u  and  states  x  that  set  the  state  derivatives  to  zero.  Initial  starting  guesses 
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Xo,  Uo  and  yo  are  given  to  the  algorithm,  where  yo  are  any  outputs  which  are  not 
also  states.  The  routine  trim  then  attempts  to  minimize  the  difference  between  the 
final  steady-state  values  and  these  initial  guesses,  while  driving  the  state  derivatives 
to  zero. 


In  this  case,  because  the  equilibrium  condition  is  steady,  level  and  unaccelerated 
flight,  the  trim  states  for  q,  p,  (j),  r,  •tp  and  /?  can  be  assumed  to  be  zero.  Level  flight 
is  assured  by  setting  the  fiight-path  angle  7  to  be  zero.  Subsequent  values  for  the 
elevator  6h,  throttle  aileron  and  rudder  Spt  and  for  the  remaining  states  V, 

a  and  6  should  exist  for  a  steady-state  solution.  The  trim  routine,  however,  does 
not  allow  the  user  to  set  any  of  the  states  or  outputs  to  a  prespecified  constant  value, 
while  then  searching  the  remaining  states,  inputs,  and  outputs  for  the  equilibrium 
values.  Instead,  trim  considers  such  values  as  “desired”  and  iterates  about  all  states 
until  equilibrium  or  some  close  proximity  to  equilibrium  is  found. 


As  a  result,  in  this  nonlinear  model  the  states  that  are  known  to  be  zero  all  take  on 
finite,  though  small  values.  The  flight-path  angle  also  is  nonzero  for  the  equilibrium 
values  selected  by  trim.  An  example  of  the  disparity  between  the  desired  and  the 
selected  values  is  shown  for  FPl  (flight  point  one)  in  Table  2.1.  The  initial  guesses 
come  from  the  trim  values  generated  by  the  Genesis  model  at  Wright-Patterson  AFB. 
Although  the  disparities  between  Genesis  and  trim  are  small  enough  to  be  sufficient 
for  the  purpose  of  deriving  a  linear  model,  the  primary  drawback  of  this  technique  is 
the  added  time  required  for  trim  to  search  all  the  states,  instead  of  just  those  which 
are  known  to  be  nonzero. 


2.4  Using  the  Nonlinear  Model 


This  section  will  give  a  brief  overview  of  the  current  capabilities  of  the  F-15  nonlinear 
simulation  model  as  well  as  the  essential  modifications  needed  for  use  with  other 
aircraft. 


Table  2.1:  Trim  Data  FPl  (9,800  ft,  0.5  M) 


Parameter  Units 


Initial  Guess 

Trim  Value 

State  Derivative 

539.08 

539.08 

-7.1636-^ 

mmm 

0.0801 

-1.8594-® 

0 

_2.79-i9 

-2.1328-5 

0.0801 

0.0803 

-2.7949-1® 

0 

-2.54-20 

2.3895-1® 

0 

-1.02-^® 

-2.8053-20 

0 

-3.22-20 

-1.8936-2® 

0 

-1.66-1® 

-3.2388-2® 

0 

-8.41-21 

-2.9103-2® 

-2.827 

-2.880 

N/A 

37.4 

37.4 

N/A 

0 

0 

N/A 

0 

0 

N/A 

0 

0.0002 

N/A 

ft/s 

rad 

rad/s 

rad 

rad/s 

rad 

rad/s 

rad 

rad 
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2.4-1  Applications 

The  nonlinear  F-15  simulation  model  is  comprised  essentially  of  the  governing  s- 
function,  which  describes  the  interconnection  of  the  key  components  which  include 
the  nonlinear  equations  of  motion,  and  the  supporting  atmospheric,  aerodynamic, 
and  propulsion  models.  Because  of  the  ease  of  implementing  an  s-function  in  the 
SIMULINK  environment,  the  model  can  be  used  for  both  the  open-  and  closed- loop 
evaluation.  The  procedure  involves  creating  an  overall  system  in  SIMULINK  in  which 
the  F-15  dynamics  are  described  by  a  single  s-function  block.  For  example,  referring 
back  to  Figure  2.3,  the  block  titled  “F-15  nonlinear  dynamics”  contains  the  open-' 
loop  s-function  f  25sf  n.  The  MATLAB  m-file  for  f  25sfn  is  found  in  Appendix  A. 


It  is  likely  that  for  other  applications,  the  user  may  desire  to  change  the  inputs, 
outputs,  or  even  the  states  of  the  original  model  set-up.  As  an  example,  such  a  change 
was  made  for  its  use  in  conjunction  with  the  TECS  control  law.  The  changes  to  these 
parameters  are  made  in  the  s-function,  as  can  be  seen  by  comparing  the  open-loop 
s-function  f 25sfn  with  the  closed-loop  TECS  s-function  f 25sfncl  in  Appendix  A. 
A  sample  of  these  changes  is  shown  below  for  clarification. 

•  Open-loop  s-function  f  25sfn  input  and  output  listings: 


ra  INPUTS  (u) 


DH  =u(l);  ra 
DD  =u(2);  m 
DA  =  u(3);  %%% 
DR  =  u(4) ;  Va 
PLAPL  =  u(5) ;  %%% 
PLAPR  =  u(6) ; 

PLASYM  =  u(7) ;  %%% 
alpdot  =  u(8) ; 


til  SYSTEM  OUTPUTS  (Y)  III 
sys(l,l)  =  x(l); 
sys(2,l)  =  x(2); 
sys(3,l)  =  x(3); 


SYMETRIC  STABILATDR  (DEG)  III 
DIFFERENTIAL  STABILATDR  (DEG)  III 
AILERON  DEFLECTION  (DEG)  III 
RUDDER  DEFLECTION  (DEG)  HI 
LEFT  PLA  (DEG)  HI 
RIGHT  PLA  (DEG)  HI 
SYMMETRIC  PLA  (DEG)  HI 
AOA  RATE  (RAD/S)  HI 

HI  V  HI 
HI  ALPHA  HI 
HI  Q  HI 


sys(4,l)  =  x(4) ; 
sys(5,l)  =  x(5); 
sys(6,l)  =  x(6) ; 
sys(7,l)  =  x(7) ; 
sys(8,l)  =  x(8) ; 
sys(9,l)  =  x(9) ; 

sys(lO,l)  =  al3+(all+al2)/(V*m*C0SBETA); 
sys(ll,l)  =  x(lO); 
sys(l2,l)  =  x(ll); 
sys(l3,l)  =  x(l2) ; 


III  THETA 
III  P 
III  PHI 
III  R 
III  PS I 
III  BETA 
III  ALPHA.DOT 
III  X-PDSITION 


III  Y-POSITION 
III  H 


•  Closed-loop  s-function  f  25sfncl  input  and  output  listings: 


III  INPUTS  (U)  111 
DH  =  u(l)+utrim(l) ; 

PLAPL  =  u(2)+utrim(2) ; 

PLAPR  =  u(3)+utrim(3) ; 

FLAPS  =  u(4)+utrim(4) ; 

DA  =  u(5)+utrim(5) ; 

DR  =  u(6)+utrim(6) ; 


III  SYMETRIC  STABILATDR  (DEG) 
III  LEFT  PLA  (DEG) 

III  RIGHT  PLA  (DEG) 

III  FLAPS  (DEG) 

III  AILERON  DEFLECTION  (DEG) 
III  RUDDER  DEFLECTION  (DEG) 


III  SYSTEM  OUTPUTS  (Y)  III 
sys(l,l)  =  (Vl+V2+V3+V4+V5)/m/g; 

sys(2,l)  =  x(l)-xx(l); 

sys(3,l)  =  x(l2)-xx(l2); 

sys(4,l)  =  x(4)-xx(4)-x(2)+xx(2); 

sys(5,l)  =  x(3)-xx(3) 

sys(6,l)  =  x(4)-xx(4) 

sys(7,l)  =  x(2)-xx(2) 

sys(8,l)  =  x(5)-xx(5) 

sys(9,l)  =  x(6)-xx(6) 

sys(lO,l)  =  be6+(bel+be2+be3+be4+be5)/(m*V) ; 

sys(ll,l)  =  x(9)-xx(9); 
sys(l2,l)  =  x(8)-xx(8); 


III  V_D0T/G 
III  V 
III  H 
III  GAMMA 

in  Q 

III  THETA 
III  ALPHA 
III  P 
III  PHI 
III  BETA_D0T 
III  BETA 
III  PSI 
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sys(l3,l)  =  r*COSPHI*SECTHETA+q*SINPHI*SECTHETA;  III  PSI_DOT  III 
sys(l4,l)  =  x(7)-xx(7);  III  R  III 

Note  that  the  closed-loop  inputs  declaration  includes  the  addition  of  the  trim  input 
initial  conditions  and  the  output  declaration  includes  the  subtraction  of  the  respec¬ 
tive  trim  initial  condition,  if  one  exists  for  that  parameter.  This  is  required  for  the 
closed-loop  s-function  since  the  values  feeding  into  the  controller  (system  outputs)  and 
the  values  of  the  controls  taken  out  of  the  controller  (system  inputs)  are  perturbations. 

If  the  user  desires  to  run  the  model  simulation  with  a  particular  set  of  initial 
conditions,  then  in  most  cases  this  is  best  done  by  declaring  these  parameters  as  Xo 
and  Uo  vectors  in  the  s-function.  This  is  particularly  appropriate  whenever  other 
blocks  in  the  overall  SIMULINK  system  contain  integrators  of  first  or  higher  order, 
such  as  in  the  case  of  a  shaping  filter.  An  excerpt  from  an  s-function  where  the  initial 
conditions  are  listed  internally  is: 

function  [sys,xO,uO]  =  f25sfncl(t ,x,u,flag) ; 
if  flag  ==  0 

III  SYSTEM  CHARACTERISTICS/INITIAL  CONDITIONS  III 
sys  =  [12  0  13  8  12  0] ; 

III  Trim  Flight  Point  2  III 

xO  =  [497.31  0.18534  0.00000  0.18944  0.00000  ... 

0.00000  0.00000  0.00000  0.00000  0  0  30000]; 
uO  =  [-6.9712  000  76.750  80.679  0  0]; 

• 

However,  in  some  instances  such  as  open-loop  simulation,  it  may  be  convenient  to  set 
these  values  in  the  MATLAB  memory  and  then  include  them  in  the  call  statement 
for  the  simulation.  An  example  MATLAB  listing  is: 

III  Trim  Flight  Point  2  III 

xO  =  [497.31  0.18534  0.00000  0.18944  0.00000  ... 
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0.00000  0.00000  0.00000  0.00000  0  0  30000]; 

IVL  Call  Simulation  tVh 
[t,x,y]  =  euler( 'f 25sim’ ,x0) ; 

Note  that  in  the  above  case,  the  initial  conditions  for  the  controls  must  be  included 
in  the  s-function  as  shown  in  the  first  example  or  as  a  value  in  a  step-input  block 
in  the  SIMULINK  system.  The  function  euler  identifies  the  method  of  numerical 
integration  the  simulation  will  use,  and  f  25sim  is  the  overall  SIMULINK  system 
within  which  the  s-function  is  a  single  block. 

2.4-2  Modification  for  Other  Aircraft 

Because  of  the  modular  design  of  this  nonlinear  simulation  model,  modification  for 
other  aircraft  is  made  relatively  easy — the  user  need  only  supply  the  modified  aero¬ 
dynamic  and  the  propulsion  models.  The  proper  format  for  these  modules  is  found 
from  the  examples  shown  in  Appendix  D.  A  summary  of  the  required  parameters  to 
be  returned  from  each  module  is  shown  below. 


•  Aerodynamic  module: 


Cl 

Coefficient  of  Lift 

EL 

Total  Rolling  Moment 

Cd 

Coefficient  of  Drag 

EM 

Total  Pitching  Moment 

Cy 

Coefficient  of  Side  Force 

EAT 

Total  Yawing  Moment 

•  Propulsion  module: 

Xx  Thrust  in  X  Body-Axis  Direction 
Yx  Thrust  in  Y  Body- Axis  Direction 

Zx  Thrust  in  Z  Body-Axis  Direction 

The  nonlinear  equations  of  motion  do  not  change  regardless  of  the  aircraft,  so  the 
actual  s-function  should  require  only  slight  modifications.  It  is  even  possible  to  use 
linearized  aerodynamic  and  propulsion  models  while  still  utilizing  the  same  basic 
equations  of  motion  for  the  nonlinear  vehicle  dynamics. 
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Many  of  the  system  parameters,  including  the  aircraft  dimensions,  the  products 
of  inertia,  and  the  moments  of  inertia,  are  defined  in  the  vector  A,  which  is  read  into 
the  aerodynamic,  propulsion,  and  atmospheric  modules  of  the  simulation.  This  vector 
is  also  used  to  store  updated  values  for  the  various  atmospheric,  aerodynamic,  and 
propulsion  parameters  as  well  as  the  states  and  control  inputs.  The  A- vector  origi¬ 
nally  came  from  the  Genesis  simulation  and  contained  2,000  parameters.  However, 
most  of  the  parameters  are  either  not  used  or  are  specific  to  the  Genesis  simulation. 
A  listing  of  the  significant  A- vector  parameters  is  shown  in  Appendix  E.  A  sample 
of  commands  for  running  both  the  open-  and  closed-loop  simulations  in  MATLAB  is 
shown  in  Appendix  F. 


Chapter  3 

MODEL  EVALUATION  AND  LINEARIZATION 


In  this  section,  the  nonlinear  F-15  model  will  be  analyzed  by  evaluating  its  open- 
loop  characteristics  in  a  steady  unaccelerated  level  flight  condition.  Aircraft  responses 
based  on  the  initial  states  derived  from  the  MATLAB  trim  function  will  be  compared 
with  those  previously  generated  by  the  Genesis  simulation  at  Wright-Patterson.  Some 
of  the  challenges  with  the  model  will  also  be  discussed,  namely  the  limitations  found 
in  the  model’s  accuracy  throughout  the  flight  envelope  and  the  difficulties  encountered 
in  determining  a  trim  point  at  a  specified  flight  condition. 

The  model  analysis  will  also  include  an  evaluation  of  the  linearized  models  at  two 
flight  points.  State  responses  to  a  pulse-input  independently  applied  to  the  elevator, 
rudder,  and  aileron  control  will  be  compared  with  those  of  the  nonlinear  model. 

3.1  Open-Loop  Nonlinear  Model  at  Trim 

Open-loop  nonlinear  simulation  is  used  to  validate  the  selected  trim  point  for  each  of 
the  flight  conditions.  Once  validated,  the  trim  points  will  be  used  to  construct  the 
linearized,  state-space  models  to  be  used  in  the  TECS  control-law  design. 

3.1.1  Comparison  to  Genesis  Simulation 

In  determining  the  trim  points  for  the  two  flight  conditions,  the  “initial  guess”  used 
for  the  aircraft  trim  states  was  obtained  from  the  Genesis  model.  Values  obtained 
from  the  MATLAB  trim  function  differ  slightly  from  those  obtained  by  the  Genesis 
simulation,  as  shown  in  Table  3.1. 

Results  of  these  variations  in  terms  of  aircraft  responses  are  shown  in  Figures 
3. 1-3. 2,  with  the  complete  set  of  plots  for  both  flight  points  shown  in  Appendix  B.  In 
the  case  of  the  first  flight  condition  (9800  ft,  0.5  M),  the  MATLAB  initial  conditions 
result  in  only  very  small  deviations  from  the  equilibrium.  The  largest  of  these  is  the 
altitude  /i,  with  an  increase  of  15  feet  over  a  time  interval  of  100  seconds,  primarily  due 
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to  a  nonzero  flight  path  angle  {a  ^  9,  Table  3.1).  Responses  from  the  Genesis  initial 
conditions  show  significant  oscillations  and  a  large  deviation  from  the  equilibrium. 
Since  the  oscillations  dampen  with  time,  this  response  is  considered  to  be  stable. 

The  second  flight  condition  (30,000  ft,  0.5  M)  shows  larger  deviations  in  the 
responses  with  the  MATLAB  initial  conditions.  The  pitch  angle  9  decreases  approxi¬ 
mately  0.1°  and  h  increases  approximately  180  feet  over  a  period  of  100  seconds.  Both 
the  velocity  V  and  the  pitch  rate  q  also  show  small  deviations  from  the  equilibrium 
values.  With  the  Genesis  initial  conditions,  larger  deviations  from  trim  are  observed, 
with  the  exception  of  the  altitude  variable.  The  Genesis  results  produce  oscillations 
with  a  significantly  higher  frequency  than  those  obtained  from  the  MATLAB  trim 
condition. 


Table  3.1:  Comparison  of  Trim  Conditions 


States/ 

FP  1 

FP  1 

FP  2 

FP  2 

Outputs 

Units 

MATLAB 

Genesis 

MATLAB 

Genesis 

V 

ft/s 

539.08 

539.08 

497.32 

497.32 

a 

rad 

0.0801 

0.0801 

0.1853 

0.1891 

9 

rad/s 

0.0 

0.0 

0.0 

0.0 

9 

rad 

0.0803 

0.0801 

0.1898 

0.1891 

P 

rad/s 

0.0 

0.0 

0.0 

0.0 

rad 

0.0 

0.0 

0.0 

0.0 

r 

rad/s 

0.0 

0.0 

0.0 

0.0 

rad 

0.0 

0.0 

0.2197 

0.0 

rad 

0.0 

0.0 

0.0 

0.0 

Sh 

deg 

-2.827 

-2.880 

-6.970 

-7.059 

^PLA 

deg 

37.39 

37.40 

78.83 

78.95 

Sa 

deg 

0.0 

0.0 

0.0 

0.0 

Sr 

deg 

0.0 

0.0 

0.0 

0.0 
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Figure  3.1:  Flight  Point  1 — Comparison  of  Aircraft  Responses  to  Initial  Conditions 
Set  at  Trim  Values:  MATLAB  (solid  line)  and  Genesis  (dashed  line). 
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Figure  3.2:  Flight  Point  2 — Comparison  of  Aircraft  Responses  to  Initial  Conditions 
Set  at  Trim  Values:  MATLAB  (solid  line)  and  Genesis  (dashed  line). 


31 


3.1.2  Limitations  in  the  Flight  Envelope 

The  two  flight  points  selected  for  the  evaluation  of  the  nonlinear  model  are  at  the 
same  Mach  number  (0.5  M).  Although  it  is  desireable  to  evaluate  the  model  over 
as  much  of  the  flight  envelope  as  possible,  the  higher  Mach  number  regime  created 
problems  for  the  nonlinear  model.  As  seen  in  Figure  3.3,  at  an  altitude  of  9,800  feet 
and  Mach  number  of  0.9  {V  —  970.34  ft/s),  high-frequency  oscillations  take  place 
after  approximately  30  seconds  of  simulation.  It  is  possible  that  these  oscillations  are 
the  results  of  inaccurate  aerodynamic  table  data.  Because  the  purpose  of  the  model 
is  to  provide  a  generic  framework  for  many  other  aircraft,  problems  specific  to  the 
F-15  nonlinear  modeling  will  be  left  for  future  investigation. 


Figure  3.3;  High-Frequency  Oscillations  at  Higher  Mach  Numbers  (9,800  ft,  0.9  M) 
3.1.3  Limitations  in  Determining  the  Trim  Point 

Although  the  challenges  associated  with  the  MATLAB  trim  function  were  already 
addressed  in  the  previous  chapter,  an  additional  problem  became  evident  when  eval¬ 
uating  the  aircraft  equilibrium  responses.  As  mentioned  in  Section  3.1.1,  the  trim 
routine  was  run  using  initial  guesses  based  on  the  trim  values  from  the  Genesis  rou¬ 
tine. 
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Clearly,  this  method  would  not  be  of  practical  use  to  other  aircraft  models,  since  such 
accurate  estimates  of  the  trim  conditions  are  usually  not  available. 

In  order  to  test  the  reliability  of  the  trim  routine,  starting  values  for  the  control 
settings  were  taken  to  be  different  from  the  original  Genesis  values,  while  keeping  the 
same  values  for  the  velocity  V,  the  angle  of  attack  a,  and  the  pitch  angle  0.  Results  of 
this  test  for  flight  point  2,  along  with  the  previous  results  obtained  using  the  Genesis 
initial  guesses,  are  shown  in  Table  3.2.  Trim  values  resulting  from  the  arbitrary  initial 
guesses  are  significantly  different  from  those  with  the  Genesis  values,  especially  for 
the  pitch  angle  0,  the  flight  path  angle  7,  and  the  throttle  setting  Spla-  Note  also 
that  trim  has  not  been  achieved  at  level  flight,  based  on  the  nonzero  values  for  7. 


Table  3.2:  Evaluation  of  Trim  Reliability  (FP2) 


Initial  Guess 

Trim 

Initial  Guess 

Trim 

Parameter 

Units 

(Genesis) 

Values 

(Arbitrary) 

Values 

R 

ft/s 

497.32 

497.32 

497.32 

a 

rad 

0.1891 

0.1853 

0.1891 

<1 

rad/s 

0.0 

0.0 

0.0 

0.0 

e 

rad 

0.1891 

0.1898 

0.1891 

0.0713 

p 

rad/s 

0.0 

0.0 

0.0 

0.0 

rad 

0.0 

0.0 

0.0 

0.0 

r 

rad/s 

0.0 

0.0 

0.0 

0.0 

rad 

0.0 

0.2197 

0.0 

0.0237 

/? 

rad 

0.0 

0.0 

0.0 

0.0 

Sh 

deg 

-7.059 

-6.970 

0.0 

-7.123 

^PLA 

deg 

78.95 

78.83 

30.00 

30.05 

Sa 

deg 

0.0 

0.0 

0.0 

0.0 

Sr 

deg 

0.0 

0.0 

0.0 

0.0 

7 

rad 

0.0 

0.0045 

0.0 

-0.1170 

Although  the  initial  guesses  for  the  elevator  of  zero  degrees  and  for  the  throttle  of 
30°  are  significantly  different  from  the  Genesis  values,  the  elevator  trim  value  of 
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Figure  3.4:  Aircraft  Responses  at  FP2  using  Trim  Conditions  Obtained  from  an 
Arbitrary  Initial  Guess  (Table  3.2). 


—7.123°  converges  roughly  to  the  Genesis  trim  value  of  —6.970°.  The  trim  throttle 
setting,  however,  remains  virtually  unchanged  from  the  initial  guess  of  30°,  resulting 
in  an  under-powered  condition.  The  trim  function  attempts  to  correct  for  this  by 
lowering  the  aircraft  nose,  allowing  the  aircraft  to  maintain  velocity  but  causing  the 
level  flight  constraint  (7  =  0)  to  be  violated.  This  situation  is  evident  from  Figure 
3.4,  showing  the  aircraft  open-loop  responses  at  this  trim  condition.  In  order  for 
the  trim  routine  to  be  useful  under  such  conditions,  it  must  provide  a  capability 
to  rigidly  set  constraints  on  the  aircraft  output  variables  such  as  constant  velocity 
and  zero  flight-path  angle.  Currently,  the  trim  routine  allows  specification  of  desired 
values  for  output,  state,  and  input  variables,  but  these  values  are  not  held  constant 
throughout  the  iteration.  In  both  examples  shown  in  Table  3.2,  the  flight  path  angle 
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7  was  “set”  to  zero  in  trim,  but  in  the  final  values  both  are  nonzero  at  the  equilibrium 
condition. 


3.2  Evaluation  of  Linearized  Model 

From  the  MATLAB  trim  data  shown  previously  in  Table  3.1,  a  linearized  state-space 
model  at  each  flight  point  is  determined  using  the  MATLAB  function  linmod.  To 
validate  these  models,  open-loop  responses  of  the  nonlinear  and  linearized  models 
are  compared  for  20-second  pulse  inputs  applied  independently  to  the  elevator,  the 
aileron,  and  the  rudder  controls.  Although  small  variations  exist,  the  linearized  model 
responses  for  both  flight  points  are  reasonably  close  to  those  of  the  nonlinear  model. 
State-space  models  of  the  full  linearized  equations  of  motion  for  both  flight  points 
are  shown  in  Appendix  C.  Stability  characteristics  of  the  open-loop  linearized  models 
are  summarized  in  Table  3.3. 


3.2.1  Longitudinal  Excitation 
Elevator  Pulse  Input 

The  elevator  pulse  input  for  both  flight  conditions  is  2°  down  for  20  seconds,  resulting 
in  an  initial  loss  in  altitude  and  an  accompanying  increase  in  airspeed.  Figures  3. 5-3.6 
show  the  response  plots  for  the  velocity  V,  the  altitude  h,  the  pitch  rate  q,  and  the 
pitch  angle  6  at  flight  points  1  and  2  respectively. 

3.2.2  Lateral- Directional  Excitation 
Aileron  Pulse  Input 

The  aileron  input  is  a  20-second  pulse  of  1°  for  flight  point  1  and  a  20-second  pulse 
of  2°  for  flight  point  2.  In  both  cases,  the  left  aileron  is  up  and  the  right  aileron 
down,  producing  a  left  turn.  The  maximum  roll  angle  reached  is  approximately  50® 
for  both  flight  conditions.  Figures  3. 7-3.8  show  the  response  plots  for  the  angle  of 
sideslip  /?,  the  yaw  rate  r,  the  roll  rate  p,  and  the  roll  angle  ^  for  flight  points  1  and 
2  respectively. 
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Table  3.3:  Open-Loop  Stability  Characteristics 


Flight  Point  1 

Mode 

Eigenvalues 

Damping 

Frequency 

(rad/s) 

Altitude 

0.0000 

1.0000 

0.0000 

Spiral 

-0.0015 

1.0000 

0.0015 

Heading 

-0.0325 

1.0000 

0.0325 

Phugoid 

-0.0055  ±  0.0805f 

0.0687 

0.0807 

Roll 

-2.1211 

1.0000 

2.1211 

Dutch  Roll 

-0.4128  ±2.5913f 

0.1573 

2.6239 

Short  Period 

-1.6407  ±  2.2486f 

0.5894 

2.7835 

Flight  Point  2 

Mode 

Eigenvalues 

Damping 

Frequency 

(rad/s) 

Altitude 

0.0000 

1.0000 

0.0000 

Spiral 

-0.0035 

1.0000 

0.0035 

Heading 

-0.0452 

1.0000 

0.0452 

Phugoid 

-0.0072  ±  0.0905f 

0.0798 

0.0908 

Roll 

-0.7086 

1.0000 

0.7086 

Dutch  Roll 

-0.3461  ±1.9067f 

0.1786 

1.9378 

Short  Period 

-0.7522  ±1.4587f 

0.4583 

1.6412 

Rudder  Pulse  Input 

A  pulse  input  to  the  twin  rudders  of  1®  right  is  applied  for  20  seconds  at  both  flight 
conditions.  The  result  is  a  right  turn,  with  a  maximum  roll  angle  for  flight  conditions 
1  and  2  of  approximately  70®  and  80°  respectively.  Figures  3.9-3.10  show  the  response 
plots  for  the  sideslip  angle  /?,  the  yaw  rate  r,  the  roll  rate  p,  and  the  roll  angle  cj)  at 
flight  points  1  and  2  respectively. 
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Figure  3.5:  Flight  Point  1 — Aircraft  Responses  to  a  20-second  Elevator  Pulse  of  2°: 
Nonlinear  Model  (solid  line)  and  Linearized  Model  (dashed  line). 
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Figure  3.6:  Flight  Point  2 — Aircraft  Responses  to  a  20-second  Elevator  Pulse  of  2°: 
Nonlinear  Model  (solid  line)  and  Linearized  Model  (dashed  line). 
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Figure  3.7:  Flight  Point  1 — Aircraft  Responses  to  a  20-second  Aileron  Pulse  of  1°: 
Nonlinear  Model  (solid  line)  and  Linearized  Model  (dashed  line). 
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Figure  3.8:  Flight  Point  2 — Aircraft  Responses  to  a  20-second  Aileron  Pulse  of  2°: 
Nonlinear  Model  (solid  line)  and  Linearized  Model  (dashed  line). 
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Figure  3.9:  Flight  Point  1 — Aircraft  Responses  to  a  20-second  Rudder  Pulse  of  1°: 
Nonlinear  Model  (solid  line)  and  Linearized  Model  (dashed  line). 
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Figure  3.10:  Flight  Point  2 — Aircraft  Responses  to  a  20-second  Rudder  Pulse  of  1®: 
Nonlinear  Model  (solid  line)  and  Linearized  Model  (dashed  line). 


Chapter  4 

LONGITUDINAL  CONTROL  USING  TECS 


4.1  Background 

Since  the  early  years  of  aircraft  automatic  flight  control  system  design,  the  funda¬ 
mental  approach  to  autopilot  design  has  evolved  around  the  single-loop  feedback 
structure.  Although  this  method  has  proven  satisfactory  for  the  time,  recently  it 
has  been  found  that  the  inherent  limitations  of  this  design  method  are  holding  back 
the  development  of  more  capable  flight  control  systems.  Extensive  studies  such  as 
discussed  in  [3]  demonstrate  that  the  conventional  design  methods  have  reached  their 
fundamental  limits.  Future  improvements  therefore  rely  on  the  development  of  new 
system  architectures. 

The  Total  Energy  Control  System  (TECS)  was  developed  in  response  to  this 
observation.  The  conventional  flight  control  system  is  fundamentally  limited  by  its 
single-loop  development  of  the  throttle-  and  elevator- command  loops,  which  neglects 
the  cross-coupling  effects  of  the  longitudinal  dynamics.  The  TECS  design,  however, 
has  the  fundamental  objective  of  integrating  flight  path  and  speed  control  based  on 
“energy  compensation”  techniques.  In  the  conventional  design,  flight-path  control  is 
achieved  exclusively  through  feedback  to  the  elevator  while  speed  control  is  achieved 
using  feedback  to  the  throttles.  One  of  the  most  serious  design  deficiencies  of  this 
method  is  the  potential  for  adverse  cross-coupling  between  the  elevator  and  throttle 
controls  once  both  feedback  loops  are  closed.  Both  stability  and  performance  of 
the  entire  system  may  be  adversely  affected.  Similar  design  deficiencies  exist  for 
conventional  lateral  control  systems. 

By  focusing  on  integration  between  the  flight  path  and  speed  control  systems, 
the  TECS  design  has  incorporated  the  following  features  to  overcome  this  limitation: 
the  use  of  a  multiloop  control  structure  with  crossfeed  paths  between  the  speed  and 
flight  path  to  both  the  throttle  and  the  elevator,  the  use  of  total  energy  rate  and 
total  energy  distribution  rate  in  the  control  design,  and  the  use  of  command  paths 
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for  control  of  flight  path  and  longitudinal  acceleration.  The  reader  is  referred  to 
the  following  sources  for  further  background  in  the  methods  of  total  energy  control 
[2,  4,  7,  8]. 

4.2  Development  of  the  TECS  Concept 

Numerous  methods  exist  for  simultaneous  operation  of  the  throttle  and  elevator  using 
total  energy  control  methods.  One  such  method  is  based  on  the  understanding  that 
the  throttle  is  fundamentally  related  to  energy  rate,  in  that  it  either  increases  or 
decreases  the  overall  system  energy.  The  elevator,  on  the  other  hand,  is  related 
to  energy  distribution.  It  causes  an  exchange  between  kinetic  and  potential  energy 
in  the  system  by  causing  the  aircraft  to  either  climb  or  descend.  An  example  of 
the  advantage  offered  by  the  TECS  method  is  with  the  simultaneous  commands 
of  increased  flight-path  angle  and  decreased  speed.  In  a  single-loop  design,  both 
the  elevator  and  throttle  would  respond,  possibly  resulting  in  over- control  or  the 
two  controls  “opposing”  one  another.  Using  the  total  energy  control  concept,  the 
system  recognizes  that  a  simultaneous  climb  and  decrease  in  airspeed  is  essentially 
an  exchange  of  the  aircraft  kinetic  energy  for  the  potential  energy,  and  would  therefore 
rely  primarily  on  the  elevator  to  achieve  the  desired  response. 

In  mathematical  terms,  the  total  energy  control  concept  is  fundamentally  based 
on  the  total  energy  E(t)  of  the  aircraft.  The  following  summarizes  the  derivation  of 
the  total  energy  control  concept  in  [9].  For  a  point  mass,  E{t)  is  described  in  terms  of 
the  aircraft  mass  m,  the  total  velocity  V{t),  the  altitude  h{t),  and  the  gravitational 
acceleration  g  as 

E{t)  — +  mgh{t)  (4.1) 

By  differentiating  equation  (4.1)  with  respect  to  time,  the  total  energy  rate  E{t)  is 

E(i)  =  msV(t)  +  7(()i  (4.2) 

where  the  flight-path  angle  7(t)  is  in  radians.  The  thrust  required  Treg(t)  is  then 
derived  from  the  equations  of  motion  along  the  flight  path.  By  assuming  the  initial 
thrust  offsets  the  drag,  and  that  variations  in  drag  are  slow,  then  the  thrust  required 
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for  a  particular  level  of  total  energy  rate  is 

Treq{t)=^^^  =  mg{^^  +  -l{t^  (4-3) 

From  equation  (4.3),  it  is  evident  that  the  aircraft  total  energy  state  is  determined 
by  the  thrust  required,  which  is  controlled  by  the  throttle.  However,  the  capacity  to 
distribute  energy  between  its  kinetic  and  potential  states  is  not  well-handled  by  the 
throttle.  Because  the  elevator  causes  relatively  little  drag  while  controlling  the  air¬ 
craft  angle  of  attack  for  small  inputs,  it  is  an  energy-conserving  mechanism.  Rather 
than  adding  to  the  total  system  energy  as  with  the  throttle,  the  elevator  essentially 
distributes  the  energy  between  the  kinetic  and  the  potential  states  while  maintain¬ 
ing  the  total  energy  constant.  Thus,  by  using  proportional  and  integral  control  on 
the  total  energy  rate  and  the  energy  distribution  rate,  the  throttle  and  the  elevator 
commands  are 

StXs)  =  rng  {Ktp  + 

where  the  flight-path  error  7e(s)  and  the  acceleration  error  ^(■s)  are 

7,(S)  = 

t{s)  =  l/(s)-Re(5) 

The  proportional  feedback  gains  to  the  throttle  and  the  elevator  are  Ktp  and  Kep, 
and  the  integral  feedback  gains  are  Kti  and  Ket  respectively. 

The  overall  closed-loop  system  using  TECS  is  shown  in  Figure  4.1.  Because  the 
model  is  designed  to  be  used  for  the  generic  aircraft,  the  inputs  have  been  reduced 
to  elevator,  throttle,  flaps,  rudder,  and  aileron.  In  the  case  of  the  F-15,  the  elevator 
is  represented  by  the  symmetric  stabilator  and  the  throttle  by  the  symmetric  PLA. 
The  differential  stabilator  and  the  individual  PLA  settings  are  left  out  to  maintain 
generality.  Flaps  are  not  defined  in  the  F-15  dynamic  model.  Inputs  to  the  longi¬ 
tudinal  TECS  controller  are  normalized  acceleration  Vc{s)/g^  velocity  R,  altitude  h, 


(4.4) 

(4.6) 
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flight-path  angle  7,  pitch  rate  q  and  pitch  angle  6.  Inputs  to  the  lateral  TECS  con¬ 
troller  are  roll  rate  p,  roll  angle  sideslip  rate  sideslip  angle  /?,  yaw  angle  0,  yaw 
rate  r  and  heading  rate  of  change  The  specific  TECS  design  for  the  longitudinal 
dynamics  will  be  discussed  next.  For  more  information  on  the  lateral  TECS  design, 
the  reader  is  referred  to  [9]. 


Figure  4.1:  Closed-Loop  Model  Block  Diagram  for  TECS  Control  Law. 


4.3  Longitudinal  TECS  Structure 

The  longitudinal  TECS  controller  is  composed  of  commands  to  altitude  he,  flight-path 
angle  7c,  normalized  acceleration  Vc/g,  and  velocity  I4.  The  controller  can  operate 
with  he  and  14  specified  only,  from  which  the  7c  and  Vc/g  values  are  derived  using 
the  proportional  gains  Kh  and  K„  respectively.  Otherwise,  the  values  for  7c  and 
Vc/g  can  be  specified  directly.  The  structure  of  the  longitudinal  TECS  controller  is 
shown  in  Figure  4.2.  Note  that  the  integral  action  on  the  airspeed  and  altitude  are 
accomplished  through  feedback  of  the  velocity  and  altitude  errors  to  the  acceleration 
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and  flight-path  commands,  respectively. 

A  state-space  representation  for  the  longitudinal  TECS  controller,  shown  in  the 
format  of  [9],  is 


I  ^  Clang 

y  ^Ciang 


^Ciang  ^ Clang  "h  ^CiangVs  long 
^^long^^long  '^^long'y  ^long 


where  the  controller  state  vector  the  measurement  vector  ysiang  the  control 

input  vector  are 


^Ciang  =  Ci  {t) ,  X c{t)f 

Vsiong  =  [A7(^),AF(<)/5,A?(t),A0(t),Ay(t),A/i(t)]^ 
^^long  -  [A5f^t),  A^ec(0]^ 
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and 


^long 


B. 


^tong 


a 


^long 


0  0 
0  0 

-1-10  0  Kh 

1-10  0  -Kh 

Kgw  Kti  0 
0  KcasKei 


D 


^long 


Kgw  Ktp  Kgw  ^tp  0  0  0  0 

KcasKep  —KcasKep  Kg  as  Kg  K gas  Kg  0  0 


The  parameters  Ktp,  Kti,  Kep  a-nd  Kei  are  the  proportional  and  integral  gains 
on  the  throttle  and  the  elevator  controls.  The  gains  Kg  and  Kg  essentially  form  the 
stability  augmentation  system  (SAS),  while  the  gains  K^  and  Kh  provide  feedback 
correction  on  the  airspeed  and  altitude  errors,  respectively.  The  gains  Kcas  and  Kgw 
are  scheduled  according  to  the  aircraft  calibrated  airspeed  Vcas  and  gross  weight  W, 
where  W  =  mg.  Determination  of  the  above  gains  will  be  accomplished  using  the 
constrained  parameter  optimization  program  SANDY.  The  SIMULINK  model  for  the 
longitudinal  TECS  controller  is  shown  in  Figure  4.3.  Implementation  of  the  controller 
with  the  F-15  nonlinear  model  is  the  topic  of  the  next  section. 


Chapter  5 

TECS  CONTROLLER  PERFORMANCE 


The  primary  objective  of  this  section  is  to  evaluate  a  control  law  designed  for 
the  F-15  using  the  linearized  aircraft  model  and  then  to  validate  this  control  law  by 
implementing  it  on  the  nonlinear  model.  In  this  example,  the  control  law  used  is  the 
longitudinal  TECS  controller  and  the  gains  have  been  selected  using  the  numerical 
optimization  program  SANDY.  Describing  the  controller  design  is  not  the  intent  of 
this  section.  The  gains  selected  for  each  flight  point  are  shown  in  Appendix  D. 

5.1  Linearized  Closed-Loop  Model  Evaluation 

Validity  of  the  linearized  models  at  the  two  flight  conditions  was  established  in  Chap¬ 
ter  3.  Now  we  will  evaluate  the  closed-loop  system  using  the  linearized  aircraft  dy¬ 
namics. 

5.1.1  Closed-Loop  Characteristics 

Stability  of  the  system  is  evaluated  initially  from  the  closed-loop  eigenvalues.  Table 

5.1  shows  the  closed- loop  eigenvalues  of  lowest  damping  and  frequency  for  each  of 
the  flight  conditions.  For  both  flight  points,  the  minimum  damping  is  0.7  and  the 
closest  pole  to  the  origin  is  at  —0.100.  A  stable  system  should  therefore  result.  The 
complete  set  of  closed- loop  eigenvalues  is  shown  in  Appendix  D. 

Robustness  of  the  system  is  evaluated  for  single-loop  stability  margins.  The  el¬ 
evator  and  throttle  control  loops  are  “broken”  independently  and  the  resulting  gain 
and  phase  margins  are  determined.  Figures  5. 1-5.2  show  the  Bode  plots  with  the 
elevator-loop  and  throttle-loop  broken  respectively  for  the  two  flight  conditions.  As 
seen  in  Table  5.2,  the  lowest  gain  and  phase  margins  substantially  exceed  even  the 
conservative  minimum  requirements  of  ±6  dB  and  ±45°  respectively,  indicating  that 
adequate  robustness  exists. 


Phase  (Degrees)  Magnitude  (dB) 
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Table  5.1:  Closed-Loop  Stability  Characteristics 


Flight  Eigenvalues 

Point 

Damping 

Frequency 

(rad/s) 

1  -0.357  ±  0.364i 

0.700 

0.510 

-0.100 

1.000 

0.100 

2  -0.124  ±0.126* 

0.700 

0.176 

-0.100 

1.000 

0.100 

Table  5.2:  Single-Loop  Stability  Margins 


Flight 

^DH  G.M. 

8dh  P.M. 

8pLA  G.M. 

SpLA  P.M 

Point 

(dB) 

(deg) 

(dB) 

(deg) 

1 

25.82 

-64.3 

oo 

-60.9 

2 

-44.19 

-62.3 

oo 

-62.3 

Frequency  (Rad/Sec)  Frequency  (Rad/Sec) 


Figure  5.1:  Elevator  Control-Loop  Bode  Plots. 
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Figure  5.2:  Throttle  Control-Loop  Bode  Plots. 


The  ability  of  the  closed-loop  system  to  reject  disturbances  is  also  evaluated. 
As  shown  in  Table  5.3,  the  root-mean-square  (rms)  values  of  the  load  factor  at  the 
aircraft  center-of-gravity  the  altitude  h,  the  total  velocity  V,  the  symmetric 

stabilator  6h  and  the  throttle  Spla  to  turbulence  of  unit  intensity  in  both  the  u-  and 
w-directions  are  all  satisfactorily  small  at  both  design  flight  points.  The  specification 
for  rizca  ^  commercial  aircraft  is  typically  a  maximum  rms  of  0.1  g.  Although  this 
requirement  is  satisfied  here,  it  would  likely  be  further  relaxed  for  a  fighter  aircraft 
such  as  the  F-15. 


Table  5.3:  RMS  Responses  to  Turbulence 


Flight 

h 

^DH 

^PLA 

Point 

(g) 

(ft) 

(ft/sec) 

(deg) 

(deg) 

1 

0.0143 

4.93 

0.0717 

0.0432 

0.2300 

2 

0.0092 

6.17 

0.0878 

0.0793 

0.7964 
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Command  bandwidths  for  the  velocity  and  altitude  variables  are  also  evaluated 
to  determine  if  adequate  control  power  exists  for  both  the  elevator  and  throttle. 
Figure  5.3  shows  the  Bode  plots  for  the  V/Vcmd  and  hjhcmd  transfer  functions  at 
flight  points  1  and  2.  The  VfVcmd  bandwidth  is  0.060  rad/s  for  flight  point  1  and 
0.058  rad/s  for  flight  point  2.  The  h/hcmd  bandwidth  is  0.060  rad/s  for  both  flight 
conditions.  Comparing  the  bandwidths  of  V/Vcmd  to  hlhcmd,  the  values  for  flight 
point  1  are  nearly  identical  while  those  for  flight  point  2  differ  by  0.002  rad/s.  In 
both  cases,  the  differences  in  bandwidth  are  small  enough  that  control  power  should 
be  equally  distributed  to  both  the  elevator  and  throttle. 

Based  on  the  above  analysis,  the  closed-loop  linearized  system  command  responses 
should  exhibit  satisfactory  stability,  robustness,  and  disturbance  rejection.  These 
responses  will  be  evaluated  in  the  next  section. 


Figure  5.3:  Vc  and  he  Command  Frequency  Responses. 
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5.1.2  Command  Responses 

The  closed-loop  command  responses  to  aircraft  velocity  and  altitude  are  determined 
using  the  SIMULINK  system  shown  in  Chapter  4,  with  the  linearized  F-15  dynamics. 
Responses  for  the  velocity  V,  altitude  h,  elevator  6fj  and  throttle  8pi^  are  discussed 
below.  The  complete  set  of  responses  are  shown  in  Appendix  D. 


Velocity  Command 

The  linearized  aircraft  responses  to  a  20  ft/s  velocity  command  are  shown  for  both 
flight  conditions  in  Figures  5.4-5.5.  The  velocity  tracking  is  satisfactory,  with  a 
settling  time  (Tg)  of  under  50  seconds  in  both  cases.  At  flight  point  1,  the  excursion 
from  the  trim  altitude  is  approximately  6  ft,  whereas  flight  point  2  is  approximately 
47  ft.  Both  the  elevator  and  throttle  control  responses  for  flight  point  2  are  more 
than  twice  those  of  flight  point  1.  Note  also  the  slight  oscillation  in  elevator  which 
occurs  in  both  conditions  at  approximately  10  seconds. 

Altitude  Command 

The  linearized  aircraft  responses  to  a  1000  ft  altitude  command  are  shown  for  both 
flight  conditions  in  Figures  5. 6-5.7.  The  altitude  tracking  shows  almost  a  10-second 
delay,  resulting  in  a  Tg  of  greater  than  60  seconds  for  both  flight  points.  However, 
the  velocity  deviation  from  trim  is  kept  very  small.  This  indicates  possibly  too  little 
authority  in  elevator  control.  Again,  flight  point  1  has  control  responses  of  less  than 
half  those  of  flight  point  2. 

Combined  Velocity  and  Altitude  Commands 

The  linearized  aircraft  responses  to  simultaneous  commands  of  20  ft/s  in  velocity 
and  1000  ft  in  altitude  are  shown  for  both  flight  conditions  in  Figures  5. 8-5.9.  The 
most  notable  difference  here  from  the  two  independent  commands  is  the  increase  in 
the  velocity  Tg  to  greater  than  50  seconds.  For  flight  point  2,  the  velocity  response 
nearly  overshoots  the  command  tracking  at  approximately  20  seconds.  The  same 
relationship  in  control  response  magnitudes  for  the  two  flight  conditions  exists  here. 
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Figure  5.4;  Flight  Point  1 — Linear  Aircraft  Responses  to  a  20  ft/s  Velocity  Command. 


Figure  5.5:  Flight  Point  2 — Linear  Aircraft  Responses  to  a  20  ft/s  Velocity  Command. 


Elevator  (deg)  Velocity  (ft/s)  I  Elevator  (deg)  Velocity  (ft/s) 


Elevator  (deg)  Velocity  (ft/s)  »  ^  Elevator  (deg) 


time  (sec) 
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;ure  5.8:  Flight  Point  1 — Linear  Aircraft  Responses  to  a  20  ft /s  Velocity  Command 


1  1000  ft  Altitude  Command. 


Figure  5.9:  Flight  Point  2 — Linear  Aircraft  Responses  to  a  20  ft/s  Velocity  Command 
and  1000  ft  Altitude  Command. 
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5.2  Nonlinear  Closed-Loop  Model  Evaluation 

Using  the  same  TECS  controller  gains  and  the  system  structure  as  Section  5.1,  the 
linearized  F-15  dynamics  are  now  replaced  by  the  nonlinear  F-15  model.  The  system 
responses  are  again  evaluated  in  the  SIMULINK  environment.  Nonlinear  system 
responses  for  the  velocity  V,  the  altitude  h,  the  elevator  6}j  and  the  throttle  Spla  a.re 
shown  below  and  compared  with  the  equivalent  linearized  responses  from  the  previous 
section.  The  complete  set  of  nonlinear  responses  is  shown  in  Appendix  D. 


5.2.1  Command  Responses 
Velocity  Command 

Nonlinear  aircraft  responses  to  a  20  ft/s  velocity  command  are  shown  for  both  flight 
conditions  in  Figures  5.10-5.11  .  Several  differences  exist  in  these  responses  from  the 
linearized  equivalents  in  Figures  5.4-5.5  .  Both  flight  points  show  velocity  responses 
with  settling  times  greater  than  50  seconds  and  with  some  overshoot.  The  altitude 
responses  exhibit  steady-state  errors  of  approximately  40  ft  and  50  ft  for  flight  points 
1  and  2  respectively.  Most  notable,  however,  are  the  high-frequency  oscillations  in 
elevator  during  the  first  40  seconds  of  the  response. 


Altitude  Command 

Nonlinear  aircraft  responses  to  a  1000  ft  altitude  command  are  shown  for  both  flight 
conditions  in  Figures  5.12-5.13  .  Here,  the  differences  between  the  nonlinear  and 
linearized  responses  are  magnified  even  further.  The  altitude  tracking  is  similar  to 
that  of  the  linear  model  until  the  command  value  is  reached — the  nonlinear  model 
exhibits  a  steady-state  error  of  approximately  50  ft  at  both  flight  points.  The  ele¬ 
vator  again  oscillates  at  high  frequency,  in  this  case  throughout  the  entire  command 
response.  The  throttle  responses  differ  significantly  from  the  linearized  responses.  A 
much  larger  steady-state  value  exists — approximately  18“  power  level  angle  for  both 
flight  conditions  compared  to  2°  and  less  than  1°  for  flight  points  1  and  2  respectively 
in  the  linear  model  responses. 
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Combined  Velocity  and  Altitude  Commands 

The  linearized  aircraft  responses  to  simultaneous  commands  of  20  ft/s  in  velocity 
and  1000  ft  in  altitude  are  shown  for  both  flight  conditions  in  Figures  5.14-5.15.  Of 
the  three  sets  of  nonlinear  responses,  the  combined  command  results  in  the  great¬ 
est  similarity  to  the  linear  model  responses.  Again,  the  elevator  response  exhibits 
high-frequency  oscillations,  but  considerably  less  than  the  altitude  command  alone. 
The  altitude  tracking  also  has  a  steady-state  error  which  is  less  than  that  of  the  alti¬ 
tude  command  alone.  The  throttle  response,  though  jerky,  corresponds  more  to  the 
linearized  throttle  response  than  either  of  the  previous  command  responses. 

5.0.1  Analysis 

The  objective  of  designing  a  satisfactory  TECS  controller  using  a  linearized  model 
derived  from  the  overall  nonlinear  system  has  been  achieved.  The  implementation  of 
the  TECS  controller  with  the  nonlinear  system,  however,  uncovered  several  shortcom¬ 
ings.  The  high-frequency  oscillations  in  the  elevator  responses  and  the  irregularities 
in  the  throttle  responses  are  most  noteworthy.  Possible  explanations  abound,  includ¬ 
ing  the  possibility  of  numerical  errors  in  the  simulation.  Time  steps  as  small  as  0.001 
seconds  were  used  in  an  effort  to  alleviate  the  high-frequency  characteristics. 
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Figure  5.10:  Flight  Point  1 — Nonlinear  Aircraft  Responses  to  a  20  ft/s  Velocity  Com 
mand. 
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Figure  5.11:  Flight  Point  2 — Nonlinear  Aircraft  Responses  to  a  20  ft/s  Velocity  Com¬ 
mand. 
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ure  5.14:  Flight  Point  1 — Nonlinear  Aircraft  Responses  to  a  20  ft/s  Velocity  Com- 


nd  and  1000  ft  Altitude  Command. 


Figure  5.15:  Flight  Point  2 — Nonlinear  Aircraft  Responses  to  a  20  ft/s  Velocity  Com 
mand  and  1000  ft  Altitude  Command. 


Chapter  6 

CONCLUSIONS 


6.1  Summary 

The  objective  of  this  report  is  to  develop  a  nonlinear  aircraft  simulation  for  the  F-15, 
use  a  subsequent  linearized  model  to  develop  a  control  law,  and  then  validate  the  final 
controller  on  the  nonlinear  model.  This  model  is  based  on  the  nonlinear,  six-degree- 
of-freedom  equations  of  motion,  with  internal  modules  representing  the  nonlinear 
aerodynamic  and  propulsion  characteristics  of  the  aircraft  and  the  atmospheric  data. 

In  validating  the  nonlinear  model  at  equilibrium  conditions,  several  shortcomings 
were  found  in  the  MATLAB  trim  function.  Nevertheless,  using  the  data  provided 
from  the  Genesis  simulation  as  an  initial  guess,  adequate  trim  conditions  were  deter¬ 
mined.  Linearization  of  the  model  at  the  two  selected  flight  conditions  was  success¬ 
fully  validated  by  the  comparison  of  the  pulse  responses  for  the  linear  and  nonlinear 
models.  Finally,  the  longitudinal  TECS  control  law  was  designed  for  the  system  using 
the  linearized  dynamics.  Although  not  ideal,  the  closed-loop  characteristics  of  the  lin¬ 
earized  model  were  certainly  adequate — satisfactory  closed-loop  stability,  robustness, 
turbulence  rejection,  and  control  bandwidths  were  all  verified.  The  linearized  com¬ 
mand  responses  exhibited  adequate  settling  times  and  good  tracking  characteristics, 
with  zero  steady-state  error. 

Implementation  of  the  TECS  controller  on  the  nonlinear  model  was  successful  with 
respect  to  achieving  stable  command  responses  similar  to  those  of  the  linearized  model 
responses.  However,  several  differences  were  identified  which  could  not  be  readily  ex¬ 
plained.  They  include  high-frequency  elevator-response  characteristics,  steady-state 
errors  in  altitude,  irregular  throttle-response  characteristics  and  high  steady-state 
throttle  values.  Since  the  intent  of  using  the  nonlinear  model  is  to  validate  the  con¬ 
troller,  these  characteristics  ideally  would  be  explained  by  shortcomings  in  the  control 
law.  However,  that  is  not  likely  the  case  here.  Based  on  the  irregular  nature  of  the 
variations  between  the  linear  and  nonlinear  model  responses,  the  greatest  likelihood 
is  that  the  interaction  between  the  nonlinear  model  and  the  controller  is  not  yet  per- 
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fected.  Future  improvements  in  the  nonlinear  model  will  be  required  to  make  it  a 
useful  tool  to  the  designer  for  the  final  validation  of  a  control-law  design. 

6.2  Recommendations  for  Future  Study 

Many  of  the  areas  requiring  further  investigation  have  already  been  mentioned.  The 
following  list  summarizes  the  primary  areas  of  concern  which  must  be  addressed  to 
make  the  nonlinear  model  a  viable  tool  to  the  controls  engineer: 

•  Improve  or  redesign  the  MATLAB  trim  function  to  allow  fixed  values  to  be  set 
for  the  states,  outputs,  or  inputs  and  to  allow  the  function  greater  flexibility  in 
terms  of  selecting  an  initial  guess.  Also,  improve  the  numerical  efficiency  of  the 
function  to  decrease  the  time  required  to  generate  a  trim  condition. 

•  Address  the  problems  with  the  high  Mach  number  flight  regime  and  verify  that 
they  are  specific  to  the  F-15  model  data  and  not  to  the  generic  nonlinear  model 
structure. 

•  Investigate  and  determine  the  causes  of  the  undesireable  closed-loop  charac¬ 
teristics  in  the  nonlinear  model  responses  mentioned  in  the  previous  section. 
Verify  whether  these  shortcomings  result  from  the  aerodynamic  and  propulsion 
modules,  making  them  specific  to  the  F-15  data,  or  if  they  are  the  result  of 
shortcomings  in  the  overall  nonlinear  model  structure. 

•  Apply  the  nonlinear  model  framework  to  other  aircraft  using  their  specific 
nonlinear  aerodynamic  and  propulsion  characteristics.  Design  an  appropriate 
control-law  using  a  linearized  model  and  then  validate  with  the  nonlinear  model. 

The  ultimate  viability  of  the  nonlinear  model  will  depend  on  the  satisfaction  of  most, 
if  not  all,  of  the  above  areas. 
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Appendix  A 

F-15  NONLINEAR  SIMULATION  S-FUNCTIONS 


A.l  S-Function  for  Open-Loop  F-15  Model 

The  following  is  a  listing  of  the  MATLAB  s-function  f 25sfn  used  for  the  open-loop 
nonlinear  F-15  simulation. 


Ill  NONLINEAR  OPEN-LOOP  F-15  SIMULATION  S-FUNCTION  HI 


function  [sys,xO]  =  f 25sfn(t,x,u,flag) 


global  A  lA 

global  V  alpha  q  theta  p  phi  r  psi  beta  xp  yp  h 
global  DH  DD  DA  DR  PLAPL  PLAPR  PLASYM  alpdot 
global  XX  utrim 


Ht  UPDATE  A,IA- 

ARRAYS  HI 

A(829) 

=  V; 

A (9 14)  =  alpha; 

A(862) 

=  q; 

A(713) 

=  h; 

A(715)  =  xp; 

A(716) 

=  yp; 

A(943) 

=  theta; 

A(86l)  =  p; 

A (942) 

=  phi; 

A(863) 

=  r; 

A(944)  =  psi; 

A(915) 

=  beta; 

A(1402) 

=  DH; 

A(1403)  =  DD; 

A(1401) 

=  DA; 

A (1404) 

=  DR; 

A(1416)  =  PLAPL; 

A(1417) 

=  PLAPR; 

A (1418)  =  PLASYM; 

HI  STANDARD  ATMOSPHERE  HI 
[AMCH.RHO.QBAR.G]  =  atmos(A); 

HI  A-ARRAY  UPDATE  HI 
A(825)  =  AMCH;  A(670)  =  RHO; 


67 


A(669)  =  QBAR;  A(772)  =  G; 

III  STABILITY  AXIS  FORCES  AMD  MOMENTS  III 

[CLFT,CD,CY,CL,CM,CN,FAX,FAY,FAZ,ALM,AMM,ANM]  =  f 25aero (A , lA) ; 

II  TRANSFER  TO  A-ARRAY 

A(1410)  =  CLFT;  A(1411)  =  CD;  A(1412)  =  CY; 

A(1413)  =  CL;  A(1414)  =  CM;  A(1415)  =  CN; 

A(748)  =  FAX;  A(749)  =  FAY;  A(750)  =  FAZ; 

A(733)  =  ALM;  A(734)  =  AMM;  A(735)  =  ANM; 

III  CALCULATE  PROPULSION  FORCES  III 

[FPX,FPY,FPZ,DCL,DCM,DCN,TAUL,TAUR,PLAL,PLAR,C0UT1C,C0UT2C, . . . 
FIRST]  =  f25eng(A,IA); 

III  TRANSFER  TO  THE  A  AND  lA  ARRAYS  III 

A(751)  =  FPX;  A(752)  =  FPY;  A(753)  =  FPZ; 

A(736)  =  DCL;  A(737)  =  DCM;  A(738)  =  DCN; 

A(1419)  =  TAUL;  A(1420)  =  TAUR; 

A(1431)  =  PLAL;  A(1432)  =  PLAR; 

A(667)  =  COUTIC;  A(668)  =  C0UT2C; 

IA(502)  =  FIRST; 

III  A, lA- ARRAY  VAR  NAMES  III 

S  =  A(659);  QBAR  =  A(669);  g  =  A(772) ; 

CLFT  =  A(1410);  CD  =A(1411);  CY  =A(1412); 

ALM  =  A(733);  AMM  =  A(734);  ANM  =  A(735) ; 

FPX  =  A(75l);  FPY  =  A(752);  FPZ  =  A(753) ; 

Ix  =  A(634);  ly  =  A(635);  Iz  =  A(636) ; 

Ixz  =  A(637);  Ixy  =A(638);  lyz  =  A(639) ; 

II  =  A(632);  12  =  A(631);  13  =  A(630) ; 

14  =  A(629);  15  =  A(628);  16  =  A(627) ; 

Dx  =  A(626);  Dy  =  A(625);  Dz  =  A(624) ; 
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detl  =  A(633);  m  =  A(658)/G; 

III  CALCULATE  UPDATED  AERO  VALUES  VH 
L  =  QBAR*S*CLFT; 

D  =  QBAR*S*CD; 

Y  =  QBAR*S*CY; 

SumL  =  ALM; 

SumM  =  AMM; 

SumN  =  ANM; 

XT  =  FPX; 

YT  =  FPY; 

ZT  =  FPZ; 

til  ANGLE  CALCULATIONS  III 
CQSTHETA  =  cos (theta); 

S INTHETA  =  sin (theta) ; 

TANTHETA  =  tan(theta) ; 

SECTHETA  =  1/cos (theta) ; 

COSBETA  =  cos (beta); 

SINBETA  =  sin(beta) ; 

TANBETA  =  tan (bet a); 

COSALPHA  =  cos (alpha); 

SINALPHA  =  sin(alpha); 

COSPHI  =  cos (phi); 

SINPHI  =  sin(phi); 

COSPSI  =  cos(psi) ; 

SINPSI  =  sin(psi) ; 

III  EQUATIONS  OF  MOTION  III 

VI  =  -D*COSBETA+Y*SINBETA+XT*COSALPHA*COSBETA; 
V2  =  YT*SINBETA+ZT*SINALPHA*COSBETA; 

V3  =  -m*g*(COSALPHA*COSBETA*S INTHETA); 

V4  =  -m>i'g*(-SINBETA*SINPHI*COSTHETA); 
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V5  =  -m*g*(-SINALPHA*CDSBETA*COSPHI*COSTHETA) : 
all  =  -L+ZT*COSALPHA-XT*S INALPHA; 

al2  =  m=t=g*(COSALPHA*CQSPHI*COSTHETA+SINALPHA*SINTHETA); 

al3  =  q-TANBETA*(p*CQSALPHA+r*SINALPHA); 

ql  =  SumL*I2+S\imM*I4+SumN*I5-p''2*(lxz*I4-Ixy*I5)  ; 

q2  =  p*q*(Ixz*I2-Iyz*I4-Dz*I5)-p*r*(Ixy*I2+Dy*I4-Iyz*I5) ; 

q3  =  q''2*(lyz*I2-Ixy*I5)-q*r*(Dx=Kl2-Ixy*I4+Ixz*I5)  ; 

q4  =  -r“2*(Iyz*I2-Ixz*I4) ; 

pi  =  SumL*Il+S'umM*I2+SumN*I3-p"2*(lxz*I2-Ixy*I3)  ; 
p2  =  p*q*(Ixz*Il-Iyz*I2-Dz*I3)-p*r*(Ixy*Il+Dy*I2-Iyz*I3) ; 
p3  =  q"2*(Iyz*Il-Ixy*I3)-q*r*(Dx*Il-Ixy*I2+Ixz*I3) ; 
p4  =  -r''2*(lyz*Il-Ixz*I2) ; 

rl  =  SumL*I3+SiMM*I5+SumN*I6-p''2*(lxz*I5-Ixy*I6)  ; 
r2  =  p*q’t=(lxz*I3-Iyz*I5-Dz*I6)-p*r*(lxy*I3+Dy*I5-Iyz*I6) ; 
r3  =  q''2*(lyz*I3-Ixy*I6)-q*r*(Dx*I3-Ixy*I5+Ixz*I6)  ; 
r4  =  -r''2*(lyz*I3-Ixz*I5) ; 

bel  =  D*SINBETA+Y*COSBETA-XT*COSALPHA>i=SINBETA; 

be2  =  YT’i'COSBETA-ZT’KSINALPHA^SINBETA; 

be3  =  m*g*(COSALPHA*SINBETA*SINTHETA); 

be4  =  m*g*(COSBETA*SINPHI*COSTHETA) ; 

be5  =  m*g*(-SINALPHA*SINBETA*COSPHI*COSTHETA); 

be6  =  p*SINALPHA-r*COSALPHA; 

xpl  =  COSALPHA*COSBETA*COSTHETA*COSPSI; 

xp2  =  SINBETA*(SINPHI*SINTHETA*CDSPSI-COSPHI*SINPSl) ; 

xp3  =  SINALPHA*COSBETA*(COSPHI*SINTHETA*SINPSI-SINPHI*COSPSl) ; 

ypl  =  COSALPHA*COSBETA*COSTHETA*SINPSI; 

yp2  =  SINBETA*(COSPHI*COSPSI+SINPHI*SINTHETA*SINPSl) : 

yp3  =  SINALPHA*COSBETA*(COSPHI*SINTHETA*SINPSI-SINPHI*COSPSl); 

hi  =  COSALPHA*COSBETA*SINTHETA; 

h2  =  SINBETA*SINPHI*COSTHETA; 

h3  =  SINALPHA*COSBETA*CQSPHI*COSTHETA; 


if  flag  ==  0 
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III  SYSTEM  CHARACTERISTICS/INITIAL  CONDITIONS  III 
sys  =  [12  0  13  8  12  0] : 
xO  =  XX ; 


III  STATES  (X)  III 

III  LONGITUDINAL  STATES  HI 


V 

=  xO(i): 

HI  TOTAL  VEHICLE  VELOCITY  (FT/S) 

III 

alpha 

=  x0(2) ; 

III  ANGLE  OF  ATTACK  (RAD) 

HI 

q 

=  x0(3) ; 

HI  PITCH  RATE  (RAD/S) 

HI 

theta 

=  xO(4); 

HI  PITCH  ANGLE  (RAD) 

HI 

III  LATERAL/DIRECTIONAL  STATES  HI 

P 

=  x0(5) ; 

HI  ROLL  RATE  (RAD/S) 

HI 

phi 

=  x0(6); 

HI  ROLL  ANGLE  (RAD) 

HI 

r 

=  x0(7) ; 

HI  YA\]  RATE  (RAD/S) 

HI 

psi 

=  x0(8) ; 

HI  YAW  ANGLE  (RAD) 

HI 

beta 

=  x0(9) ; 

HI  SIDESLIP  ANGLE  (RAD) 

HI 

HI  EARTH-RELATIVE 

POSITION  STATES  HI 

xp 

=  xO(lO) 

HI  X-DIRECTION  POSITION  (FT) 

HI 

yp 

=  xO(ll) 

HI  Y-DIRECTION  POSITION  (FT) 

HI 

h 

=  x0(12) 

HI  ALTITUDE  (FT) 

HI 

elseif  abs(flag)  ==  1 
III  STATES  (X)  III 


HI  LONGITUDINAL 

STATES  HI 

V 

=  x(l) 

HI  TOTAL  VEHICLE  VELOCITY 

(FT/S)  HI 

alpha 

=  x(2) 

HI  ANGLE  OF  ATTACK  (RAD) 

HI 

q 

=  x(3) 

HI  PITCH  RATE  (RAD/S) 

HI 

theta 

=  x(4) 

HI  PITCH  ANGLE  (RAD) 

HI 

HI  LATERAL/DIRECTIONAL  STATES  HI 

P 

=  x(5) 

HI  ROLL  RATE  (RAD/S) 

HI 

phi 

=  x(6) 

HI  ROLL  ANGLE  (RAD) 

HI 

r 

=  x(7) 

HI  YAW  RATE  (RAD/S) 

HI 

psi 

=  x(8) 

HI  YAW  ANGLE  (RAD) 

HI 
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beta 

=  x(9);  III  SIDESLIP  ANGLE  (RAD) 

III 

III  EARTH- 

-RELATIVE  POSITION  STATES  III 

xp 

=  x(lO);  m  X- 

■DIRECTION  POSITION  (FT) 

III 

yp 

=  x(ll);  m  Y- 

-DIRECTION  POSITION  (FT) 

III 

h 

=  x(l2);  IXl  ALTITUDE  (FT) 

III 

III  INPUTS  (U)  III 

DH 

=  u(l)-^utrim(l)  ; 

III  SYMETRIC  STABILATOR 

(DEG) 

III 

DD 

=  u(2) ; 

III  DIFFERENTIAL  STABILATOR  (DEG) 

III 

DA 

=  u(3)-^utrim(5)  ; 

III  AILERON  DEFLECTION 

(DEG) 

III 

DR 

=  u(4)■^utriIn(6) ; 

III  RUDDER  DEFLECTION  (DEG) 

III 

PLAPL 

=  u(5)+utrim(2) ; 

III  LEFT  PLA  (DEG) 

III 

PLAPR 

=  u(6)+utrim(3) ; 

III  RIGHT  PLA  (DEG) 

III 

PLASYM 

=  u(7) ; 

III  SYMMETRIC  PLA  (DEG) 

III 

alpdot 

=  u(8) ; 

III  AOA  RATE  (RAD/S) 

III 

III  STATE 

DERIVATIVES  (dX/dT)  III 

sysCl.l)  = 

(Vl+V2+V3+V4-i-V5)/itt; 

sys(2,l)  = 

=  al3+(all-^al2)/(V*m*C0SBETA); 

sys(3,l)  =  (ql+q2+q3+q4)/detl ; 
sys(4,l)  =  q*COSPHI-r*SINPHI; 
sys(5,l)  =  (pl+p2+p3+p4)/detl; 

sys(6,l)  =  p+q*SINPHI*TANTHETA+r*COSPHI*TANTHETA; 

sys(7,l)  =  (rl+r2+r3+r4)/detl; 

sys(8,l)  =  r*COSPHI*SECTHETA+q*SINPHI*SECTHETA; 

sys(9,l)  =  be6+(bel+be2+be3+be4+be5)/(m*V) ; 

sys(lO,l)  =  V*(xpl+xp2+xp3) ; 

sys(ll,l)  =  V*(ypl+yp2+yp3); 

sys(l2,l)  =  V*(hl-li2-h3); 


elseif  flag  ==  3 

III  SYSTEM  OUTPUTS  (Y)  III 
sys(l,l)  =  x(l); 
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sys(2,l)  =  x(2) ; 
sys(3,l)  =  x(3) ; 
sys(4,l)  =  x(4) ; 
sys(5,l)  =  x(5) ; 
sys(6,l)  =  x(6) ; 
sys(7,l)  =  x(7) ; 
sys(8,l)  =  x(8) ; 
sys(9,l)  =  x(9) ; 

sys(10,l)  =  al3+(all+al2)/(V*m*C0SBETA); 
sys(ll,l)  =  x(10) ; 
sys(12,l)  =  x(ll); 
sys(13,l)  =  x(12) ; 

else 

TLl  ALL  OTHER  FLAGS  UNDECLARED  IM 
sys  =  []  ; 
end 
end 

A. 2  S-Function  for  Closed-Loop  F-15  Model 

The  following  is  a  listing  of  the  MATLAB  s-function  f  25sfncl  used  for  the  closed- 
loop  nonlinear  F-15  simulation: 

*/•/,'/  NONLINEAR  CLOSED-LOOP  F-15  SIMULATION  S-FUNCTION  ’///'/. 
function  [sys,xO]  =  f25sfncl(t,x,u,flag) 
global  A  lA 

global  V  alpha  q  theta  p  phi  r  psi  beta  xp  yp  h 
global  DH  PLAPL  PLAPR  FLAPS  DA  DR  PLASYM  DD 
global  XX  utrim 


V/X  UPDATE  A,IA-ARRAYS  %ll 


A(829) 

=  V; 

A (9 14)  =  alpha; 

A(862)  =  q; 

A(713) 

=  h; 

A(715)  =  xp; 

A(716)  =  yp; 

A(943) 

=  theta; 

A(86l)  =  p; 

A(942)  =  phi; 

A(863) 

=  r; 

A(944)  =  psi; 

A(915)  =  beta; 

A(1402) 

=  DH; 

A(1403)  =  DD; 

A(1401)  =  DA; 

A(1404) 

=  DR; 

A(1416)  =  PLAPL; 

A(1417)  =  PLAPR; 

A(1418) 

=  PLASYM; 

lit  STANDARD  ATMOSPHERE  III 
[AMCH,RHQ,QBAR,G]  =  atmos(A); 

III  A-ARRAY  UPDATE  VH 
A(825)  =  AMCH;  A(670)  =  RHQ ; 

A(669)  =  QBAR;  A(772)  =  G; 

III  STABILITY  AXIS  FORCES  AND  MOMENTS  III 

[CLFT,CD,CY,CL,CM,CN,FAX,FAY,FAZ,ALM,AMM,ANM]  =  f25aero(A,IA) ; 


II  TRANSFER  TO  THE  A-ARRAY 
A(1410)  =  CLFT;  A(1411)  = 
A(1413)  =  CL;  A(1414)  = 
A(748)  =  FAX;  A(749)  = 
A(733)  =  ALM;  A(734)  = 


in 

CD; 

A(1412)  =  CY; 

CM; 

A(1415)  =  CN; 

FAY; 

A(750)  =  FAZ; 

AMM; 

A(735)  =  ANM; 

III  CALCULATE  PROPULSION  FORCES  III 

CFPX,FPY,FPZ,DCL,DCM,DCN,TAUL,TAUR,PLAL,PLAR,C0UT1C,C0UT2C, . . . 
FIRST]  =  f25eng(A,IA); 


II  TRANSFER  TO  THE  A  AND  lA  ARRAYS  III 


A(75l)  =  FPX; 

A(736)  =  DCL; 

A(1419)  =  TAUL; 
A(1431)  =  PLAL; 


A (752)  =  FPY; 

A (737)  =  DCM; 

A(1420)  =  TAUR; 
A (1432)  =  PLAR; 


A(753)  =  FPZ; 
A(738)  =  DCN; 
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A(667)  =  COUTIC;  A(668)  =  C0UT2C; 

IA(502)  =  FIRST; 


III 

A. 

lA-ARRAY  VAR  NAMES 

III 

s 

= 

A(659);  QBAR 

= 

A(669) ; 

S 

= 

A(772); 

CLFT 

= 

A(1410);  CD 

= 

A(1411) ; 

CY 

= 

A(1412) ; 

ALM 

= 

A(733) 

AMM 

= 

A (734) ; 

ANM 

= 

A(735) 

FPX 

= 

A(751) 

FPY 

= 

A(752) ; 

FPZ 

= 

A(753) 

Ix 

= 

A (634) 

ly 

= 

A(635) ; 

Iz 

= 

A(636) 

Ixz 

= 

A(637) 

Ixy 

= 

A(638) ; 

lyz 

= 

A(639) 

11 

= 

A(632) 

12 

= 

A(631) ; 

13 

= 

A(630) 

14 

= 

A(629) 

15 

= 

A(628) ; 

16 

= 

A(627) 

Dx 

= 

A(626) 

Dy 

= 

A(625) ; 

Dz 

= 

A(624) 

detl 

= 

A(633) 

m 

= 

A(658)/G; 

til  CALCULATE  UPDATED  AERO  VALUES  III 
L  =  QBAR*S*CLFT; 

D  =  QBAR=kS*CD; 

Y  =  qBAR*S*CY; 

SumL  =  ALM; 

SumM  =  AMM; 

SumM  =  ANM; 

XT  =  FPX; 

YT  =  FPY; 

ZT  =  FPZ; 

III  ANGLE  CALCULATIONS  III 
COSTHETA  =  cos (theta); 

SINTHETA  =  sin(theta); 

TANTHETA  =  tan(theta); 

SECTHETA  =  1/cos (theta) ; 

COSBETA  =  cos (beta); 

SINBETA  =  sin(beta); 


TANBETA 

COSALPHA 

S INALPHA 

COSPHI 

SINPHI 

CQSPSI 

SIMPSI 


tan (beta) ; 
cos (alpha) ; 
sin (alpha) ; 
cos (phi) ; 
sin(phi) ; 
cos(psi) ; 
sin(psi) ; 


VH  EQUATIONS  OF  MOTION  WL 

VI  =  -D*COSBETA+Y*SINBETA+XT*COSALPHA*COSBETA; 

V2  =  YT*SINBETA+ZT*SINALPHA*COSBETA; 

V3  =  -m*g*(COSALPHA*COSBETA*SINTHETA) ; 

V4  =  -in*g*(-SINBETA*SINPHI*COSTHETA)  ; 

V5  =  -m*g*(-SINALPHA*COSBETA*COSPHI*COSTHETA); 
all  =  -L+ZT*COSALPHA-XT*S INALPHA; 

al2  =  m*g* (COSALPHA*COSPHI*COSTHETA+SINALPHA*S INTHETA) ; 

al3  =  q-TANBETA*(p*COSALPHA+r*SINALPHA); 

ql  =  SumL*I2+SumM*I4+SumN*I5-p“2*(Ixz*I4-Ixy*I5) ; 

q2  =  p*q*(Ixz*I2-Iyz*I4-Dz*I5)-p*r*(Ixy*I2+Dy*I4-Iyz*I5) ; 

q3  =  q''2*(lyz*I2-Ixy*I5)-q*r=*=(Dx*I2-Ixy*I4+Ixz*I5)  ; 

q4  =  -r''2*(Iyz*I2-Ixz*I4)  ; 

pi  =  SumL*Il+S'uraM*I2+SumN*I3-p''2*(Ixz*I2-Ixy*I3)  ; 
p2  =  p*q*(Ixz*Il-Iyz*I2-Dz*I3)-p*r*(lxy*Il+Dy*I2-Iyz*I3) ; 
p3  =  q''2*(lyz*Il-Ixy*I3)-q*r*(Dx*Il-Ixy*I2+Ixz*l3) ; 
p4  =  -r''2*(lyz=i'Il-Ixz*I2) ; 

rl  =  SnmL*I3+S'umM*I5+SumN*I6-p''2*(lxz*I5-Ixy*I6)  ; 
r2  =  p*q*(lxz*I3-Iyz*I5-Dz*I6)-p*r=i'(lxy*I3+Dy*I5-Iyz*I6) ; 
r3  =  q''2*(lyz*I3-Ixy*I6)-q*r*(Dx*I3-Ixy*I5+Ixz*I6) ; 
r4  =  -r''2*(lyz*I3-Ixz*I5) ; 

bel  =  D*SINBETA+Y*COSBETA-XT*COSALPHA*SINBETA: 
be2  =  YT*COSBETA-ZT*SINALPHA*SINBETA; 
be3  =  m*g*(COSALPHA*SINBETA*SINTHETA); 
be4  =  m*g*(COSBETA*SINPHI*COSTHETA); 
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be5  =  m*g*(-SINALPHA*SINBETA*COSPHI*COSTHETA); 

be6  =  p*SINALPHA-r*COSALPHA; 

xpl  =  COSALPHA*COSBETA*COSTHETA*COSPSI; 

xp2  =  SINBETA*(SIMPHI*SINTHETA*COSPSI-CDSPHI*SINPSl) ; 

xp3  =  SINALPHA*COSBETA*(CQSPHI*SINTHETA*SIMPSI-SIMPHI*COSPSl) ; 

ypl  =  COSALPHA*COSBETA*CQSTHETA*SINPSI; 

yp2  =  SINBETA*(COSPHI*CQSPSI+SINPHI*SINTHETA*SINPSI); 

yp3  =  SINALPHA*COSBETA*(COSPHI*SINTHETA*SINPSI-SINPHI*COSPSI) ; 

hi  =  COSALPHA*COSBETA*SINTHETA; 

h2  =  SINBETA*SINPHI*COSTHETA: 

h3  =  SINALPHA*COSBETA*COSPHI*COSTHETA; 

if  flag  ==  0 

lit  SYSTEM  CHARACTERISTICS/INITIAL  CONDITIONS  III 
sys  =  [12  0  14  6  12  0] ; 
xO  =  XX ; 


III  STATES  (X)  III 

III  LONGITUDINAL  STATES  HI 

V  =  xO(l);  III  TOTAL  VEHICLE  VELOCITY  (FT/S)  III 
alpha  =  x0(2);  III  ANGLE  OF  ATTACK  (RAD)  HI 
q  =  x0(3);  PITCH  RATE  (RAD/S)  HI 
theta  =  xO(4);  HI  PITCH  ANGLE  (RAD)  HI 
III  LATERAL/DIRECTIONAL  STATES  HI 

p  =  x0(5);  HI  ROLL  RATE  (RAD/S)  HI 
phi  =  x0(6);  HI  ROLL  ANGLE  (RAD)  HI 
r  =  x0(7);  HI  YAW  RATE  (RAD/S)  HI 
psi  =  x0(8);  y,y.y.  YAW  ANGLE  (RAD)  HI 
beta  =  x0(9);  HI  SIDESLIP  ANGLE  (RAD)  HI 
HI  EARTH-RELATIVE  POSITION  STATES  HI 

xp  =  xO(lO);  y.y,y.  X-DIRECTION  POSITION  (FT)  HI 
yp  =  xO(ll);  III  Y-DIRECTION  POSITION  (FT)  HI 
h  =  x0(l2);  III  ALTITUDE  (FT)  HI 
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III  LONGITUDINAL 

STATES  HI 

V 

=  x(l) 

HI  TOTAL  VEHICLE  VELOCITY  (FT/S) 

III 

alpha 

=  x(2) 

HI  ANGLE  OF  ATTACK  (RAD) 

III 

q 

=  x(3) 

HI  PITCH  RATE  (RAD/S) 

III 

theta 

=  x(4) 

HI  PITCH  ANGLE  (RAD) 

HI 

III  LATERAL/DIRECTIONAL  STATES  HI 

P 

=  x(5) 

HI  ROLL  RATE  (RAD/S) 

III 

phi 

=  x(6) 

HI  ROLL  ANGLE  (RAD) 

HI 

r 

=  x(7) 

HI  YAW  RATE  (RAD/S) 

III 

psi 

=  x(8) 

HI  YAV  ANGLE  (RAD) 

III 

beta 

=  x(9) 

HI  SIDESLIP  ANGLE  (RAD) 

III 

•///'/  EARTH-RELATIVE  POSITION  STATES  VLt 

xp  =  x(lO);  m  X-DIRECTION  POSITION  (FT)  IVh 

yp  =x(ll);  y,y,y,  Y-DIRECTION  POSITION  (FT)  IVh 

h  =x(l2);  y,y;/o  ALTITUDE  (FT)  VH 


III  INPUTS  (U)  III 

DH  =  u(l)+utrim(l) ;  III  SYMETRIC  STABILATOR  (DEG)  III 

PLAPL  =  u(2)+'atrim(2) ;  LEFT  PLA(DEG)  III 

PLAPR  =  u(3)+utrim(3);  '/X/,  RIGHT  PLA  (DEG)  III 

FLAPS  =  u(4)+utrim(4);  III  FLAPS  (DEG)  III 

DA  =  u(5)+utrim(5) ;  '/X/.  AILERON  DEFLECTION  (DEG)  III 

DR  =  u(6)+utriia(6) ;  RUDDER  DEFLECTION  (DEG)  III 

PLASYM  =  0 ; 

DD  =0; 

III  STATE  DERIVATIVES  (dX/dT)  III 
sys(l,l)  =  (Vl+V2+V3+V4+V5)/ni; 
sys(2,l)  =  al3+(all+al2)/(V*m*CQSBETA); 
sys(3,l)  =  (ql+q2+q3+q4)/detl ; 
sys(4,l)  =  q*COSPHI-r*SINPHI; 
sys(5,l)  =  (pl+p2+p3+p4)/detl; 
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sys (6 , 1)  =  p+q*SINPHI*TANTHETA+r*COSPHI*TANTHETA ; 

sys(7,l)  =  (rl+r2+r3+r4)/detl ; 

sys (8,1)  =  r*COSPHI*SECTHETA+q*SINPHI*SECTHETA; 

sys (9,1)  =  be6+(bel+be2+be3+be4+be5)/(m*V) ; 

sys(lO,l)  =  V*(xpl+xp2+xp3) ; 

sys (11,1)  =  V*(ypl+yp2+yp3) ; 

sys(12,l)  =  V*(hl-h2-h3) ; 

elseif  flag  ==  3 

til  SYSTEM  OUTPUTS  (Y)  III 
sys(l,l)  =  (Vl+V2+V3+V4+V5)/m/g; 

sys(2,l)  =  x(l)-xx(l); 

sys(3,l)  =  x(l2)-xx(l2)  ; 

sys(4,l)  =  x(4)-xx(4)-x(2)+xx(2) ; 

sys(5,l)  =  x(3)-xx(3); 

sys(6,l)  =  x(4)-xx(4); 

sys(7,l)  =  x(2)-xx(2); 

sys(8,l)  =  x(5)-xx(5); 

sys(9,l)  =  x(6)-xx(6); 

sys(lO,l)  =  be6+(bel+be2+be3+be4+be5)/(m*V) ; 
sys(ll,l)  =  x(9)-xx(9); 
sys(12,l)  =  x(8)-xx(8); 

sys (13,1)  =  r*CQSPHI*SECTHETA+q*SINPHI*SECTHETA; 
sys(l4,l)  =  x(7)-xx(7); 

else 

III  ALL  OTHER  FLAGS  UNDECLARED  III 
sys  =  []  ; 
end 


end 


Appendix  B 

NONLINEAR  MODEL  RESPONSES  AT  EQUILIBRIUM 

Aircraft  responses  at  the  equilibrium  point  obtained  from  the  Genesis  and  MAT- 
LAB  trim  functions  are  shown  in  Figures  B.1-B.4  for  both  flight  conditions. 
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Figure  B.l:  Flight  Point  1 — Aircraft  Responses  to  Initial  Conditions  Set  at  Trim 
Values:  MATLAB  (solid  line)  and  Genesis  (dashed  line). 
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Figure  B.2:  Flight  Point  1 — Aircraft  Responses  to  Initial  Conditions  Set  at  Trim 
Values:  MATLAB  (solid  line)  and  Genesis  (dashed  line). 
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Appendix  C 

LINEARIZED  STATE-SPACE  MODELS 


The  following  are  state-space  models  of  the  linearized  F-15  model  for  the  two 
flight  conditions  listed  in  Table  1.1. 


States  =  [Ay(/t/sec),  Aa{rad),Aq(rad/sec),  A0{rad),  Ap{rad/sec),  A^[rad), 
Ar(rad/sec),  A^(rad),  A(S(rad) ,  Ah{f  t)]^ 


Inputs  =  [ASnideg],  A6pLAL{deg),  ASpLAnideg),  ASpiapideg),  A6A{deg),A6R{deg)f 


Outputs  =  [AF /g,  AV{ftjsec),  Ah{ft),  A'y{rad),  Aq{radlsec),  A9{rad),  Aa{rad), 
Ap[radfsec),  A^(rad),  A$[radfsec),  A/?(rad),  Axl>{rad),  A^[radf  sec),  Ar(rad/sec)j 
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Flight  Point  1 

h  =  9,800  ft,  Vtas  =  539  ft/s 


-0.0137 

-3.461 

0 

-32.144 

0 

-1.28-® 

-2.17“^ 

-0.779 

1 

-2.99-® 

0 

-2.97-^ 

4.07“^ 

-5.802 

-2.501 

-6.50-^ 

3.15-® 

4.04-® 

0 

0 

1 

0 

0 

3.23-2® 

7.45-22 

2.42-1® 

-1.53-1® 

0 

-2.21 

0 

0 

0 

-8.24-2® 

-3.25-2® 

1 

-2.25-2® 

-9.44-23 

9.45-2® 

00 

00 

1 

to 

o 

0 

-7.43-2 

0 

0 

0 

-1.02-1® 

-2.60-21 

0 

-2.80-1® 

1.16-22 

-2.69-2® 

0 

4.40-21 

8.00-2 

0.0594 

2.46-^ 

-539.1 

0 

539.1 

0 

2.15-^ 

0 

0 

-1.04-® 

-8.74-®  1 

0 

0 

1.66-12 

CO 

00 

1 

00 

-3.15-® 

0 

-4.55-1® 

-3.78-® 

1.02-1® 

0 

0 

0 

1.397 

0 

-27.05 

0 

0.0805 

0 

0 

0 

-0.5739 

0 

4.675 

0 

1.003 

0 

0 

0 

-0.997 

0 

-0.192 

-1.36-2^ 

0 

0 

-6.63-^ 

0 
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-.0718 

0.1067 

0.1014 

0 

0 

0 

-1.57-3 

-1.59-® 

-1.51-® 

0 

0 

0 

-0.1508 

4.34-® 

4.13-® 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.169 

0.0266 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2.24-3 

-0.0487 

0 

0 

0 

0 

0 

0 

-1.12-24 

1.67-24 

1.58-24 

0 

-3.83-^ 

6.32-4 

0 

0 

0 

0 

0 

0 

-4.25"^ 

-0.1077 

0 

-1 

0 

-3.99-^ 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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0 

0 

-1 

0 

1 
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0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

0 
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0 
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0 

0 

0 

0 

0 

1 

1.15-22 

-2.69-2® 

0 

4.40-2^ 

0.0800 

0.0594 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1.03-^® 

-2.61-24 

0 

-2.80-4® 

0 

0 

0 

0 

0 

0 

Flight  Point  2 

/i-  30,000  ft,  7tv1s  =  497  ft/s 


-0.0200 

-38.596 

0 

-2.61-^ 

-0.355 

1 

2.65-^ 

-2.301 

-1.148 

0 

0 

1 

1.37-21 

-3.24-1® 

-9.36-2® 

0 

0 

O 

1 

o 

1 

-1.41-22 

1.74-1® 

5.34-2® 

0 

0 

-5.42-2® 

1.10-23 

-7.77-2® 

0 

4.10-3 

-497.3 

0 

0 

0 

-4.70-^ 

0 

0 

5.37-12 

-3.15-® 

0 

-4.94-1® 

5.33-2® 

0 

0 

1.210 

0 

-13.163 

0.1917 

0 

0 

-0.301 

0 

1.455 

1.018 

0 

0 

-0.983 

0 

-0.0944 

0 

0 

-1.02-® 

-32.083 

0 

-2.90-® 

-2.81-^ 

0 

-3.12-^ 

-6.40-^ 

3.15-® 

1.13-® 

0 

0 

O 

1 

00 

0 

-1.050 

0 

-2.05-2® 

1 

2.63-23 

0 

-7.52-® 

0 

-3.86-21 

0 

1.40-22 

1 

o 

1 

0.1842 

0.0634 

497.31 

0 

4.50-^ 

-2.19-^ 

8.35-® 

-7.68-^ 

0 

0 

0 

0 

0 

-1.20-2® 

0 
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-0.142 

0.0389 

0.0370 

0 

0 

0 

-7.04"^ 

-1.47"'^ 

-1.39-® 

0 

0 

0 

-0.0554 
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1.28-® 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.0559 

5.47-® 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1.40-^ 

-0.0214 

0 

0 

0 

0 

0 

0 

-7  77-24 

2.12-24 

2.02-24 

0 

2.99-® 

3.17-^ 

0 

0 

0 

0 

0 

0 

-6.24-4 

-1.203 

0 

-1 

0 

-9.05-^ 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 

1 

0 

0 

0 
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1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

1.10-2® 

-7.77-2® 

0 

-1.10-21 

0.184 

0.0633 

0 

0 
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Appendix  D 

CLOSED-LOOP  MODEL  ANALYSIS 


D.l  TECS  Controller  Gains 

The  following  is  a  summary  of  the  gains  obtained  for  the  longitudinal  TECS  controller 
by  the  optimization  routine  SANDY.  The  linearized  state-space  F-15  model  for  each 
flight  point  was  used  in  determining  the  gains.  Refer  to  Chapter  4  for  an  explanation 
of  each  gain. 


Parameter 

FlightPointl 

FlightPoint2 

Kep 

21.57 

-93.52 

Kei 

22.34 

72.73 

Ktp 

-102.7 

-159.7 

Kti 

185.0 

291.8 

lU 

-0.0094 

-0.0118 

Kh 

-1.36"^ 

-1.07-4 

Ke 

94.4008 

1000 

K, 

34.1758 

246.5 

Kgw 

1 

1 

Kcas 

1 

1 

D.2  Closed-Loop  Eigenvalues 

A  listing  of  the  closed-loop  eigenvalues  and  their  corresponding  damping  and  frequen¬ 
cies  is  given  below  for  both  flight  points. 
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Eigenvalues 

Flight  Point  1 

Damping 

Frequency(rad/s) 

-0.0998 

1.000 

0.998 

-0.0800  ±  o.oeooz 

0.800 

0.100 

-0.0800  ±  o.oeooi 

0.800 

0.100 

-0.2926  ±  0.0962i 

0.950 

0.308 

-0.3081 

1.000 

0.308 

-0.3567  ±  0.3639i 

0.700 

0.510 

-0.3672  ±  0.3748i 

0.700 

0.525 

-0.7687 

1.000 

0.769 

-5.0127  ±5.1017^■ 

0.701 

7.152 

-6.8946 

1.000 

6.895 

Eigenvalues 

Flight  Point  2 

Damping 

Frequency(rad /s) 

-0.1000 

1.000 

0.100 

-0.0800  ±  0.0600i 

0.800 

0.100 

-0.0800  ±  0.0600j 

0.800 

0.100 

-0.1235  ±0.1260i 

0.700 

0.176 

-0.2700  ±  0.0887i 

0.950 

0.284 

-0.2842 

1.000 

0.284 

-0.4308  ±  0.4395i 

0.700 

0.615 

-0.7769 

1.000 

0.777 

-4.6791 

1.000 

4.679 

-8.7309  ±  7.801h' 

0.746 

11.71 

D.3  Linearized  Model  Responses 

Figures  D.1-D.6  are  the  command  response  plots  for  the  linearized  F-15  model. 


Figure  D.l:  Flight  Point  1 — Linear  Aircraft  Responses  to  a  20  ft/s  Velocity  Com- 
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Figure  D.3:  Flight  Point  1 — Linear  Aircraft  Responses  to  a  1000  ft  Altitude  Com- 
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Figure  D.4:  Flight  Point  2 — Linear  Aircraft  Responses  to  a  1000  ft  Altitude  Com 
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Figure  D.6:  Flight  Point  2 — Linear  Aircraft 
and  1000  ft  Altitude  Command. 
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DA  Nonlinear  Model  Responses 

Figures  D.7-D.12  are  the  command  response  plots  for  the  nonlinear  F-15  model. 
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Figure  D.7:  Flight  Point  1 — Nonlinear  Aircraft  Responses  to  a  20  ft/s  Velocity  Com 
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Figure  D.IO:  Flight  Point  2 — Nonlinear  Aircraft  Responses  to  a  1000  ft  Altitude 
Command. 
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Figure  D.ll:  Flight  Point  1 — Nonlinear  Aircraft  Responses  to  a  20  ft/s  Velocity 
Command  and  1000  ft  Altitude  Command. 
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Figure  D.12:  Flight  Point  2 — Nonlinear  Aircraft  Responses  to  a  20  ft/s  Velocity 
Command  and  1000  ft  Altitude  Command, 


Appendix  E 

F-15  NONLINEAR  SIMULATION  MODULES 


E.l  A  -  Vector  Listing 

The  following  is  a  summary  of  the  important  A- vector  parameters  used  in  the  F-15 
simulation. 


A(1418) 

=  PLASYM 

VH  SYMMETRIC  PLA  (DEG) 

A(1417) 

=  PLAPR 

%tl  RIGHT  PLA  (DEG) 

•/.'/.'/. 

A(1416) 

=  PLAPL 

•/.‘/.y.  LEFT  PLA  (DEG) 

A(1415) 

=  CN 

y.y.y.  yawing  moment  coefficient 

A(1414) 

=  CM 

VH  PITCHING  MOMENT  COEFFICIENT 

A(1413) 

=  CL 

y.y.'/.  ROLLING  MOMENT  COEFFICIENT 

A(1412) 

=  CY 

y,y,y.  side  force  coefficient 

A(1411) 

=  CD 

%%t  DRAG  COEFFICIENT 

A(1410) 

=  CLFT 

W.  LIFT  COEFFICIENT 

A (1404) 

=  DR 

W,  RUDDER  DEFLECTION  (DEG) 

A(1403) 

=  DD 

y.y.y.  differential  stabilator  (deg) 

A(1402) 

=  DH 

'/.y.y.  symmetric  stabilator  (deg) 

A(1401) 

=  DA 

AILERON  DEFLECTION  (DEG) 

A(944) 

=  psi 

HEADING  ANGLE  (RAD) 

A(943) 

=  theta 

•/.'/.•/.  PITCH  ANGLE  (RAD) 

A(942) 

=  phi 

ROLL  ANGLE  (RAD) 

•/.'/.•/. 

A(917) 

=  betadot 

SIDESLIP  RATE  (RAD/S) 

til 

A(916) 

=  alp dot 

ANGLE  OF  ATTACK  RATE  (RAD/S) 

III 

A(915) 

=  beta 

'/.'/.'/.  SIDESLIP  ANGLE  (RAD) 

III 

A(914) 

=  alpha 

'/.'/.'/.  ANGLE  OF  ATTACK  (RAD) 

III 

A(863) 

=  r 

'/.'/,'/.  YAW  RATE  (RAD/S) 

III 

A(862) 

=  q 

•/.'/.'/.  PITCH  RATE  (RAD/S) 

III 

A(86l) 

=  P 

'/.'/.'/.  ROLL  RATE  (RAD/S) 

III 

A(830) 

= 

Vdot 

in 

ACCELERATION  (FT/S''2) 

A(829) 

= 

V 

in 

VELOCITY  (FT/S) 

A(825) 

= 

AMCH 

III 

MACH  NUMBER 

A(772) 

= 

G 

III 

ACCL  OF  GRAVITY  (FT/S''2) 

A(753) 

= 

FPZ 

III 

PROPULSION  FORCE  Z-AXIS  (LBS) 

A(752) 

= 

FPY 

III 

PROPULSION  FORCE  Y-AXIS  (LBS) 

A(751) 

= 

FPX 

III 

PROPULSION  FORCE  X-AXIS  (LBS) 

A(750) 

= 

FAZ 

III 

AERO  FORCE  Z-AXIS  (LBS) 

A(749) 

= 

FAY 

III 

AERO  FORCE  Y-AXIS  (LBS) 

A (748) 

= 

FAX 

III 

AERO  FORCE  X-AXIS  (LBS) 

A(735) 

= 

ANM 

in 

YAWING  MOMENT  (FT-LBS) 

A(734) 

= 

AMM 

III 

PITCHING  MOMENT  (FT-LBS) 

A(733) 

= 

ALM 

III 

ROLLING  MOMENT  (FT-LBS) 

A(716) 

= 

yp 

III 

Y-POSITION  (FT) 

A(715) 

= 

xp 

III 

X-POSITION  (FT) 

A(713) 

= 

h 

III 

ALTITUDE  (FT) 

A(670) 

= 

RHO 

III 

DENSITY  (SLUGS/FT‘3) 

A(669) 

= 

QBAR 

III 

DYNAMIC  PRESSURE  (SL/FT-S~2) 

A(661) 

cbar 

III 

WING  CHORD  (FT) 

A(660) 

= 

b 

III 

WING  SPAN  (FT) 

A(659) 

= 

S 

III 

WING  AREA  (FT* 2) 

A(658) 

= 

W 

III 

AIRCRAFT  WEIGHT  (LBS) 

A(639) 

= 

lyz 

III 

PRODUCT  OF  INERTIA  Y-Z  PLANE 

A(638) 

= 

Ixy 

III 

PRODUCT  OF  INERTIA  X-Y  PLANE 

A(637) 

= 

Ixz 

III 

PRODUCT  OF  INERTIA  X-Z  PLANE 

A(636) 

= 

Iz 

III 

MOMENT  OF  INERTIA  Z-AXIS 

A(635) 

= 

ly 

ni 

MOMENT  OF  INERTIA  Y-AXIS 

A(634) 

= 

Ix 

III 

MOMENT  OF  INERTIA  X-AXIS 
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E.2  F-15  Nonlinear  Aerodynamic  Model  Listing 

The  following  is  a  listing  of  the  modules  which  comprise  the  nonlinear  aerodynamic 
model  for  the  F-15. 


•  MATLAB  function:  f  25aero 


function  [CLFT,CD,CY,CL,CM,CN,FAX,FAY,FAZ,ALM,AMM,ANM]  =  ... 
f 25aero(A,IA) ; 

tit  DYNAMICS  (TRIM,  HOLD,  RUM,  AND  LINEARIZATION)  tit 
[CLFT,CD,CY,CL,CM,CN,FAX,FAY,FAZ,ALM,AMM,ANM]  =  ccalc(A); 


•  MATLAB  function:  ccalc 


ttt  ROUTINE  TO  COMPUTE  AERODYNAMIC  FORCE  AND  MOMENT  ttt 
ttt  COEFFICIENTS  ttt 

ttt  WRITTEN  JUNE  16,1978 
ttt  LEE  DUKE  NASA/DFRC 

ttt  MODIFIED  OCTOBER,  1990  (CLEAN-UP,  ADD  VAR  DECL) 

ttt  RANDY  BRUMBAUGH  PRC 

ttt  BUG  FIXES  12/21/90  RWB 

ttt  CONVERSION  TO  MATLAB  12/7/93  JPD 

function  [CLFT,CD,CY,CL,CM,CN,FAX,FAY,FAZ,ALM,AMM,ANM]  =  ccalc(A) 


ttt  ASSIGN  A- ARRAY  VAR  NAMES  (INDEX  +  1001)  ttt 


DGR  =  A(981);  GO 

=  A(994); 

ALP 

=  A(914); 

BTA  =  A(915) 

B 

=  A(660); 

V 

=  A(829); 

CBAR  =  A(661) 

H 

=  A(713); 

VDOT 

=  A(830); 

ALPDOT  =  A(916) 

BTADOT 

=  A(917); 

Q 

=  A(862) ; 

P  =  A(86l) 

R 

=  A(863); 

G 

=  A(772) ; 

AZ  =  A(765) 

DA 

=  A(1401); 

DH 

=  A(1402) 
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DT  =  A(1403); 
PHI  =  A(942); 

S  =  A(659); 
DELZ  =  A(655): 
AMCH  =  A(825); 
DSB45  =0.0; 


DR  =  A (1404); 

THA  =  A(943) ; 

DELX  =  A(654); 

PLAPL  =  A (1416); 


DSB  =  A(1405); 

QBAR  =  A(669); 

DELY  =A(656); 
PLAPR  =A(1417); 


tit  PLA  IS  AVERAGE  OF  ENGINE  THROTTLE  SETTINGS  ttt 
PLA  =  (PLAPL+PLAPR)/2; 

ALPD  =  ALP*DGR; 

BTAD  =  BTA*DGR; 

if  V  ~=  0 

B2V  =  B/(2.0*V); 

CB2V  =  CBAR/(2.0*V) ; 
end 

COSALP  =  cos (ALP); 

SINALP  =  sin(ALP); 

CQSTHA  =  cos (THA); 

COSPHI  =  cos (PHI); 

ALPl  =  ALPD; 

if  ALPl  <  0,  ALPl  =  0.0;  end 
if  ALPl  >  10,  ALPl  =  10.0;  end 
ALP2  =  ALPD; 

if  ALPD  <  10,  ALP2  =  10.0;  end 


LATEST  =  0;  ttt  LOGICAL  'FALSE'  ttt 
if  ALPD  >=  25 

if  abs(BTAD)  <=  10 

if  AMCH  <=  0.6,  LATEST  =1;  end  ttt  LOGICAL  'TRUE'  ttt 
end; 
end; 


DCDALT  =  0. 0005/10000.0* (30000. 0-H) ; 


no 


CDNDZK  =1.0; 

if  PLA  >  83,  CDNOZK  =  0.0;  end 


lit  DO  TABLE  LOOK-UP  III 
[FMCOEF]  =  tlu(A); 


CLFTB 

= 

FMCOEF (1) ; 

DCLNZ 

= 

FMC0EF(2);  DCLO 

= 

FMCOEF (3) ; 

DCLAl 

= 

FMCOEF (4) ; 

DCLA2 

= 

FMCOEF (5);  CDB 

= 

FMCOEF (6) ; 

DCDNOZ 

= 

FMCOEF (7) ; 

DCDSB 

= 

FMCOEF (8);  CYBl 

= 

FMCOEF (9) ; 

CYB2 

= 

FMCOEF (10); 

CYDA 

= 

FMCOEF(ll);  CYDD 

= 

FMCOEF (12) 

DCYDR 

= 

FMCOEF (13) ; 

DRYK 

= 

FMCOEF (14);  CLBl 

= 

FMCOEF (15) 

CLB2 

= 

FMCOEF (16) ; 

CLP 

= 

FMC0EF(17);  CLR 

= 

FMCOEF (18) 

CLDA 

= 

FMCOEF (19); 

CLDD 

= 

FMCOEF (20);  DCLDR 

FMC0EF(21) 

DRLK 

= 

FMCOEF (22); 

DCLSB 

= 

FMC0EF(23);  CMB 

= 

FMCOEF (24) 

DCMNZ 

= 

FMCOEF (25) ; 

DCMO 

= 

FMCOEF (26);  DNOSB 

= 

FMCOEF (27) 

DNO 

= 

FMCOEF (28); 

CMQ 

= 

FMC0EF(29);  CMAD 

= 

FMCOEF (30) 

CNBl 

= 

FMCOEF (31); 

CNB2 

= 

FMCOEF (32);  CNP 

= 

FMCOEF (33) 

CNR 

= 

FMCOEF (34) ; 

CNDA 

= 

FMC0EF(35);  CNDD 

= 

FMCOEF (36) 

CNDR 

= 

FMCOEF (37); 

DRNK 

= 

FMC0EF(38);  DCNSB 

FMCOEF (39) ; 

III  DETERMINE  TERMS  TO  BE  USED  IN  FORCE  AND  MOMENT  EQUATIONS  HI 
CYB  =  CYBl; 

if  LATEST  ==  1,  CYB  =  CYB2;  end 
CLB  =  CLBl; 
if  LATEST  ==  1 

if  PLA  <  35,  CLB  =  CLB2;  end 
end; 

CNB  =  CNBl; 

if  LATEST  ==  1,  CNB  =  CNB2;  end 
ANZC  =  (AZ-G*C0STHA*C0SPHI)/G0; 

III  COMPUTE  AERODYNAMIC  FORCE  AND  MOMENT  COEFFICIENTS  HI 
HI  RWB  ADDED  0.95  HI 

CLFT  =  0.95*CLFTB+DCLNZ*ANZC+DSB45*(DCL0+DCLA1*ALP1+DCLA2*. . . 


(ALP2-10.0)); 

CY  =  CYB+CYDA*DA+CYDD*DT-DCYDR*DRYK; 

CL  =  CLB+CLDA*DA+CLDD*DT-DCLDR*DRLK+B2V*(CLP*P+CLR*R)+. 

DCLSB*DSB45; 

CDX  =  CDB; 

if  ALPD  >  32,  CDX  =  CLFTB*SINALP/COSALP ;  end 
if  ALPD  <  40 

if  ALPD  >  32,  CDX  =  (CDX-CDB)/8.0*(ALPD-32.0)+CDB;  end 
end; 

Vhl  RWB  ADDED  1.02  */'/*/. 

CD  =  1.02*CDX+DCDALT+DCDNQZ*CDN0ZK+DSB45*DCDSB; 

til  CALCULATES  ALPHADQT  FOR  PITCHING  MOMENT  III 
CM  =  CMB+DCMNZ*ANZC+CB2V=t=(CMQ=i'Q+CMAD*ALPD0T)+CLFTB*DN0+.  . 
DSB45*(DCM0+DN0SB*CLFTB) ; 

CN  =  CNB+CNDA*DA+CNDD*DT+CNDR*DR*DRNK+B2V*(CNP*P+CNR*R)+. 
DCNSB*DSB45 ; 

III  ADDED  TO  CORRECT  FORCE  AND  MOMENT  COEFFICIENTS  FROM  III 
III  REFERENCE  POINT  TO  THE  A/C  CG.  THESE  LINES  OF  CODE  III 
III  WERE  MERGED  INTO  THIS  SUBROUTINE  FROM  THE  SUBROUTINE  III 
III  'CGCALC'  WHICH  HAS  BEEN  DELETED.  Ill 

III  COMPUTE  BODY  AXIS  FORCE  COEFFICIENTS  III 
FX  =  -CD*COSALP+CLFT*SINALP; 

FY  =  CY; 

FZ  =  -CD*SINALP-CLFT*COSALP; 

III  COMPUTE  MOMENTS  INDUCED  BY  CHANGE  OF  REFERENCE  III 
DELTL  =  (FZ*DELY-FY*DELZ)/B; 

DELTM  =  (FX>KDELZ-FZ*DELX)/CBAR; 

DELTN  =  (FY*DELX-FX*DELY)/B; 
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III  CORRECT  MOMENT  COEFFICIENTS  III 
CL  =  CL+DELTL; 

CM  =  CM+DELTM; 

CM  =  CN+DELTM; 

III  CALCULATE  FORCE  TERMS  III 
FAX  =  FX*QBAR*S: 

FAY  =  FY*QBAR*S; 

FAZ  =  FZ*QBAR*S; 

III  CALCULATE  MOMENT  TERMS  III 
ALM  =  QBAR*S*B*CL; 

AMM  =  QBAR*S*CBAR*CM; 

ANM  =  QBAR*S*B*CN; 

•  MATLAB  function:  tlu 


function  [FMCOEF]  =  tlu(A) 


III  ROUTINE  TO  DO  TABLE  LOOK-UP  FOR  AEROMODEL  III 

III  WRITTEN  JUNE  1,  1978  III 
III  LEE  DUKE  NASA/DFRC  III 
III  MODIFIED  OCTOBER  1990  III 
III  RANDY  BRUMBAUGH  PRC  III 

III  DECLARE  GLOBAL  VARIABLES  III 

global  FIOIA  F105A  F108A  F107A  F106A  F201A  F203A  F206A  F301A 
global  F302A  F313A  F314A  F315A  F316A  F401A  F402A  F411A  F412A 
global  F413A  F414A  F415A  F416A  F418A  F501A  F505A  F506A  F507A 
global  F509A  F513A  F514A  F601A  F602A  F611A  F612A  F613A  F614A 
global  F615A  F617A  F618A 


III  ASSIGN  A-ARRAY  VAR  NAMES  III 


DR  =  A (1404); 

ALP  =  A(914); 

BTA  =  A(915); 

DGR  =  A(98l); 

AMCH  =  A(825); 

DH  =  A(1402); 

III  CREATE  DATA  ARRAYS  HI 

AMCHA  =  [0.2  0.4  0.6  0.8  0.9  1.0  1.1  1.2  1.4  1.6]; 

AMCHB  =  [0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6]; 

ALPHAA  =  [-12.0  -8.0  -4.0  0.0  4.0  8.0  12.0  16.0]; 

ALPHAA  =  [ALPHAA,  20.0  24.0  28.0  32.0  36.0  40.0 
44.0  48.0  52.0  56.0  60.0]  ; 

ALPHAB  =  [25.0  30.0  35.0  40.0  45.0  50.0  55.0  60.0]; 
DHX  =  [-25.0  -15.0  -5.0  5.0  15.0]; 

CLX  =  [-1.0  -0.8  -0.6  -0.4  -0.2  0.0  0.2  0.4  0.6 
0.8  1.0  1.2  1.4  1.6] ; 

BETAA  =  [0.0  4.0  8.0  12.0  16.0  20.0  24.0  28.0]; 
BETAB  =  [-10.0  -5.0  0.0  5.0  10. O]; 

BETAC  =  [0.0  4.0  8.0  12.0  16.0]; 

BETAD  =  [10.0  20.0  30.0]  ; 

DRXl  =  [0.0  10.0  20.0  30.0]; 

ADR  =  abs(DR); 

HI  CONVERT  TO  DEGREES  HI 
ALPD  =  ALP*DGR; 

BTAD  =  BTA*DGR; 

ABSBTA  =  abs(BTAD); 

HX  CALCULATE  INDICES  INTO  ARRAYS  HI 
IDXMAl  =  1; 
for  I  =  2:9, 

if  AMCH  >=  AMCHA(I),  IDXMAl  =  I;  end 
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end; 

IDXMA2  =  IDXMAl+1; 

XDXMBl  =  AMCH/0.2+0.01: 

IDXMBl  =  floor(XDXMBl) ; 
if  IDXMBl  <  1,  IDXMBl  =  1;  end 
if  IDXMBl  >  7.  IDXMBl  =7;  end 
IDXMCl  =  1; 

XDXAAl  =  (ALPD+12.0)/4.0+1.01; 
IDXAAl  =  floor(XDXAAl); 
if  IDXAAl  <  1,  IDXAAl  =  1;  end 
if  IDXAAl  >  18,  IDXAAl  =  18;  end 
XDXABl  =  (ALPD-25.0)/5.0+1.01; 
IDXABl  =  floor(XDXABl); 
if  IDXABl  <  1,  IDXABl  =  1 ;  end 
if  IDXABl  >  7,  IDXABl  =7;  end 
XDXDHl  =  (DH+25.0)/10.0+1.01; 
IDXDHl  =  f loor(XDXDHl) ; 
if  IDXDHl  <  1,  IDXDHl  =  1;  end 
if  IDXDHl  >  4,  IDXDHl  =4;  end 
XDXBAl  =  ABSBTA/4.0+1.01; 

IDXBAl  =  floor (XDXBAl); 
if  IDXBAl  <  1,  IDXBAl  =  1;  end 
if  IDXBAl  >  7,  IDXBAl  =7;  end 
XDXBBl  =  (BTAD+10.0)/5.0+1.01; 
IDXBBl  =  floor(XDXBBl); 


if 

IDXBBl 

< 

1, 

IDXBBl 

=  1; 

end 

if 

IDXBBl 

> 

4, 

IDXBBl 

=  4; 

end 

IDXBCl  =  : 

IDXBAl 

•  f 

if 

IDXBCl 

< 

1, 

IDXBCl 

=  1; 

end 

if 

IDXBCl 

> 

4, 

IDXBCl 

=  4; 

end 

IDXBDl  =  1; 

if  ABSBTA  >  10.0,  IDXBDl  =2;  end 
XDXDRl  =  ADR/10.0+1.01; 


IDXDRl  =  floor(XDXDRl) ; 

if  IDXDRl  <  1,  IDXDRl  =  1;  end 

if  IDXDRl  >  3,  IDXDRl  =3;  end 

IMDXAl  =  IDXMA1+10*(IDXAA1-1+19*(IDXDH1-1)); 

INDXA2  =  IMDXAl+l; 

IMDXA3  =  INDXAl+10; 

INDXA4  =  INDXA2+10; 

INDXA5  =  INDXAl+190; 

INDXA6  =  INDXA2+190; 

INDXA7  =  INDXA3+190; 

IMDXA8  =  INDXA4+190; 

INDXCl  =  IDXMA1+10*(IDXAA1-1); 

INDXC2  =  INDXCl+1; 

INDXC3  =  INDXCl+10; 

IMDXC4  =  INDXC2+10; 

INDXDl  =  IDXMA1+10*(IDXAA1-1+19*(IDXBA1-1)) ; 
INDXD2  =  INDXDl+l; 

INDXD3  =  INDXDl+lO; 

INDXD4  =  INDXD2+10; 

INDXD5  =  INDXDl+190; 

INDXD6  =  INDXD2+190; 

INDXD7  =  INDXD3+190; 

INDXD8  =  INDXD4+190; 

INDXEl  =  IDXMC1+2*(IDXAB1-1+8*(IDXBB1-1)) ; 
INDXE2  =  INDXEl+1; 

INDXE3  =  INDXEl+2; 

INDXE4  =  INDXE2+2; 

INDXE5  =  INDXEl+16; 

INDXE6  =  INDXE2+16; 

INDXE7  =  INDXE3+16; 

INDXE8  =  INDXE4+16; 

INDXFl  =  IDXMA1+10*(IDXAA1-1+19*(IDXDR1-1)) ; 
INDXF2  =  INDXFl+1; 


116 


INDXF3  =  INDXFl+10; 

INDXF4  =  INDXF2+10; 

INDXF5  =  INDXFl+190; 

INDXF6  =  IMDXF2+190; 

INDXF7  =  IMDXF3+190; 

INDXF8  =  INDXF4+190; 

INDXGl  =  IDXMA1+10*(IDXAA1-1+19*(IDXBC1-1)); 

INDXG2  =  INDXGl+1; 

INDXG3  =  IMDXGl+10; 

INDXG4  =  INDXG2+10; 

INDXG5  =  INDXGl+190; 

INDXG6  =  INDXG2+190; 

INDXG7  =  INDXG3+190; 

INDXG8  =  INDXG4+190; 

INDXHl  =  IDXMA1+10*(IDXAA1-1+19*(IDXBD1-1)) ; 

INDXH2  =  INDXHl+1; 

INDXH3  =  INDXHl+10; 

INDXH4  =  INDXH2+10; 

INDXH5  =  INDXHl+190; 

INDXH6  =  INDXH2+190; 

INDXH7  =  INDXH3+190; 

INDXH8  =  INDXH4+190; 

lit  COMPUTE  INTERPOLATION  RATIOS  III 

RATOMA  =  (AMCH-AMCHA(IDXMA1))/(AMCHA(IDXMA2)-AMCHA(IDXMA1)) ; 
RATOMB  =  (AMCH-AMCHB(lDXMBl))/0.2; 

RATOMC  =  (AMCH-0.2)/0.4; 

RATOAA  =  (ALPD-ALPHAA(lDXAAl))/4.0; 

RATOAB  =  (ALPD-ALPHAB(lDXABl))/5.0; 

RATODH  =  (DH*0.8-DHX(IDXDH1))/10.0; 

RATOBA  =  (ABSBTA-BETAA(lDXBAl))/4.0; 

RATOBB  =  (BTAD-BETAB(IDXBBl))/5.0; 

RATOBC  =  (ABSBTA-BETAC(IDXBCl))/4.0; 
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RATOBD  =  (ABSBTA-BETAD(IDXBD1))/10.0; 
if  ABSBTA  <  10.0,  RATOBD  =0.0;  end 
RATODR  =  (ADR-DRX1(IDXDR1))/10.0; 


if 

AMCH 

> 

1.6, 

RATOMA 

= 

1.0; 

end 

if 

AMCH 

> 

1.6, 

RATOMB 

= 

1.0; 

end 

if 

AMCH 

< 

0.2, 

RATOMC 

= 

0.0; 

end 

if 

AMCH 

> 

0.6, 

RATOMC 

= 

1.0; 

end 

if 

ALPD 

< 

-12.0, 

RATOAA 

= 

0.0; 

end 

if 

ALPD 

> 

60.0, 

RATOAA 

= 

1.0; 

end 

if 

ALPD 

< 

25.0, 

RATOAB 

= 

0.0; 

end 

if 

ALPD 

> 

60.0, 

RATOAB 

= 

1.0; 

end 

if 

ABSBTA 

> 

28.0, 

RATOBA 

= 

1.0; 

end 

if 

BTAD 

< 

-10.0, 

RATOBB 

= 

0.0; 

end 

if 

BTAD 

> 

10.0, 

RATOBB 

= 

1.0; 

end 

if 

ABSBTA 

> 

16.0, 

RATOBC 

= 

1.0; 

end 

if 

ABSBTA 

> 

30.0, 

RATOBD 

= 

1.0; 

end 

RAIMA 1  =  l.O-RATOMA; 

RATMCl  =  l.O-RATOMC; 

RATAAl  =  1.0-RATQAA; 

RATABl  =  l.O-RATOAB; 

RATDHl  =  l.O-RATODH; 

RATBAl  =  l.O-RATOBA; 

RATBBl  =  l.O-RATOBB; 

RATBCl  =  l.O-RATOBC; 

RATBDl  =  l.O-RATOBD; 

RATDRl  =  1.0 -RATODR; 

IVh  TABLE  LOOK  UP  III 

FlOlAl  =  RAT0MA*F101A(INDXA2)+F101A(INDXA1)*RATMA1 
F101A2  =  RAT0MA*F101A(INDXA4)+F101A(INDXA3)*RATMA1 
F101A3  =  RAT0MA*F101A(INDXA6)+F101A(INDXA5)*RATMA1 
F101A4  =  RAT0MA*F101A(INDXA8)+F101A(INDXA7)*RATMA1 
FIOIBI  =  RAT0AA*F101A2+F101A1*RATAA1; 
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F101B2  =  RAT0AA*F101A4+F101A3*RATAA1; 

FIOIBI  =  RATDAA*F101A2+F101A1*RATAA1; 

CLFTB  =  RATDDH*F101B2+F101B1*RATDH1; 

DCLNZ  =  RATOMA*F105A(IDXMA2)+F105A(IDXMA1)*RATMA1; 
DCLO  =  RAT0MA*F106A(IDXMA2)+F106A(IDXMA1)*RATMA1; 
DCLAl  =  RATQMA*F107A(IDXMA2)+F107A(lDXMAl)=t<RATMAl; 
DCLA2  =  RATQMA*F108A(IDXMA2)+F108A(IDXMA1)*RATMA1; 
DCDNOZ  =  RATQMA*F203A(IDXMA2)+F203A(IDXMA1)*RATMA1; 
DCDSB  =  RATQMA*F206A(IDXMA2)+F206A(IDXMA1)*RATMA1; 
F301A1  =  RAT0MA*F301A(INDXD2)+F301A(INDXD1)*RATMA1; 
F301A2  =  RATQMA*F301A(INDXD4)+F301A(INDXD3)*RATMA1; 
F301A3  =  RAT0MA*F301A(INDXD6)+F301A(INDXD5)*RATMA1; 
F301A4  =  RAT0MA*F301A(INDXD8)+F301A(INDXD7)*RATMA1; 
F301B1  =  RAT0AA*F301A2+F301A1*RATAA1; 

F301B2  =  RAT0AA*F301A4+F301A3*RATAA1; 

F301  =  RAT0BA*F301B2+F301B1*RATBA1; 

CYBl  =  F301*sign(BTA) ; 

F302A1  =  RAT0MC*F302A(INDXE2)+F302A(INDXE1)*RATMC1; 
F302A2  =  RAT0MC*F302A(INDXE4)+F302A(INDXE3)*RATMC1; 
F302A3  =  RAT0MC*F302A(INDXE6)+F302A(INDXE5)*RATMC1; 
F302A4  =  RAT0MC*F302A(INDXE8)+F302A(INDXE7)*RATMC1; 
F302B1  =  RAT0AB*F302A2+F302A1*RATAB1; 

F302B2  =  RAT0AB*F302A4+F302A3*RATAB1; 

CYB2  =  RAT0BB*F302B2+F302B1*RATBB1; 

F313A1  =  RAT0MA*F313A(INDXC2)+F313A(INDXC1)*RATMA1; 
F313A2  =  RAT0MA*F313A(INDXC4)+F313A(INDXC3)*RATMA1; 
CYDA  =  RAT0AA*F313A2+F313A1*RATAA1; 

F314A1  =  RAT0MA*F314A(INDXC2)+F314A(INDXC1)*RATMA1; 
F314A2  =  RAT0MA*F314A(INDXC4)+F314A(INDXC3)*RATMA1; 
CYDD  =  RAT0AA*F314A2+F314A1*RATAA1; 

F315A1  =  RAT0MA*F315A(lNDXF2)+F315A(INDXFl)*RA’fflAl: 
F315A2  =  RAT0MA*F315A(INDXF4)+F315A(INDXF3)*RA™A1; 
F315A3  =  RAT0MA*F315A(INDXF6)+F315A(INDXF5)*RATMA1; 


F315A4  =  RAT0MA*F315A(INDXF8)+F315A(INDXF7)*RATMA1; 
F315B1  =  RAT0AA*F315A2+F315A1*RATAA1; 

F315B2  =  RAT0AA*F315A4+F315A3*RATAA1; 

F315  =  RAT0DR*F315B2+F315B1*RATDR1; 

DCYDR  =  F315*sign(DR) ; 

DRYK  =  RAT0MA*F316A(IDXMA2)+F316A(IDXMA1)*RATMA1; 
F401A1  =  RAT0MA*F401A(INDXD2)+F401A  (INDXDl)=t=RATMAl 
F401A2  =  RAT0MA*F401A(INDXD4)+F401A  (INDXD3)*RATMA1 
F401A3  =  RAT0MA*F401A(INDXD6)+F401A  (INDXD5)*RATMA1 
F401A4  =  RAT0MA*F401A(INDXD8)+F401A  (INDXDr)*RATMAl 
F401B1  =  RATQAA*F401A2+F401A1*RATAA1; 

F401B2  =  RAT0AA*F401A4+F401A3*RATAA1; 

F401  =  RAT0BA*F401B2+F401B1*RATBA1; 

CLBl  =  F401*sign(BTA); 

F402A1  =  RATQMC*F402A(INDXE2)+F402A  (INDXE1)*RATMC1 
F402A2  =  RATQMC*F402A(INDXE4)+F402A  (INDXE3)*RATMC1 
F402A3  =  RAT0MC*F402A(INDXE6)+F402A  (INDXE5)*RATMC1 
F402A4  =  RAT0MC*F402A(INDXE8)+F402A  (INDXE7)*RATMC1 
F402B1  =  RAT0AB*F402A2+F402A1*RATAB1; 

F402B2  =  RATOAB+F402A4+F402A3+RATAB1; 

CLB2  =  RAT0BB*F402B2+F402B1*RATBB1; 

F411A1  =  RAT0MA*F411A(INDXC2)+F411A(INDXC1)*RATMA1; 
F411A2  =  RAT0MA*F411A(INDXC4)+F411A(INDXC3)*RATMA1; 
CLP  =  RAT0AA*F411A2+F411A1*RATAA1; 

F412A1  =  RAT0MA*F412A(INDXC2)+F412A(INDXC1)*RATMA1; 
F412A2  =  RATQMA*F412A(INDXC4)+F412A(INDXC3)*RATMA1; 
CLR  =  RAT0AA*F412A2+F412A1*RATAA1; 

F413A1  =  RAT0MA*F413A(INDXC2)+F413A(INDXC1)*RATMA1; 
F413A2  =  RAT0MA*F413A(INDXC4)+F413A(INDXC3)*RATMA1; 
CLDA  =  RAT0AA*F413A2+F413A1*RATAA1; 

F414A1  =  RAT0MA*F414A(INDXC2)+F414A(INDXC1)*RATMA1: 
F414A2  =  RAT0MA*F414A(INDXC4)+F414A(INDXC3)*RATMA1; 
CLDD  =  RAT0AA*F414A2+F414A1*RATAA1; 
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F415A1  =  RAT0MA=t=F415A(lMDXF2)+F415A(lNDXFl)=t=RATMAl; 
F415A2  =  RAT0MA*F415A(INDXF4)+F415A(INDXF3)*RATMA1; 
F415A3  =  RAT0MA*F415A(INDXF6)+F415A(INDXF5)*RATMA1; 
F415A4  =  RAT0MA*F415A(INDXF8)+F415A(INDXF7)*RATMA1; 
F415B1  =  RAT0AA*F415A2+F415A1*RATAA1; 

F415B2  =  RAT0AA*F415A4+F415A3*RATAA1; 

F415  =  RAT0DR*F415B2+F415B1*RATDR1; 

DCLDR  =  F415*sign(DR); 

DRLK  =  RATDMA*F416A(IDXMA2)+F416A(IDXMA1)*RATMA1; 
F418A1  =  RAT0MA*F418A(INDXG2)+F418A(INDXG1)*RATMA1; 
F418A2  =  RATQMA*F418A(INDXG4)+F418A(INDXG3)*RATMA1; 
F418A3  =  RATQMA*F418A(IMDXG6)+F418A(INDXG5)*RATMA1; 
F418A4  =  RAT0MA*F418A(INDXG8)+F418A(INDXG7)*RATMA1; 
F418B1  =  RATOAA=t<F418A2+F418Al*RATAAl; 

F418B2  =  RAT0AA*F418A4+F418A3*RATAA1; 

DCLSB  =  RAT0BC*F418B2+F418B1*RATBC1; 

F501A1  =  RAT0MA*F501A(INDXA2)+F501A(INDXA1)*RATMA1; 
F501A2  =  RAT0MA*F501A(INDXA4)+F501A(INDXA3)*RATMA1; 
F501A3  =  RAT0MA*F501A(INDXA6)+F501A(INDXA5)*RATMA1; 
F501A4  =  RAT0MA*F501A(INDXA8)+F501A(INDXA7)*RATMA1; 
F501B1  =  RAT0AA*F501A2+F501A1*RATAA1; 

F501B2  =  RAT0AA*F501A4+F501A3*RATAA1; 

F501B1  =  RAT0AA*F501A2+F501A1*RATAA1; 

CMB  =  RAT0DH*F501B2+F501B1*RATDH1; 

DCMNZ  =  RAT0MA*F505A(IDXMA2)+F505A(IDXMA1)*RATMA1; 
DCMO  =  RAT0MA*F506A(lDXMA2)+F506A(lDXMAl)>t=RATMAl; 
DNOSB  =  RAT0MA*F507A(IDXMA2)+F507A(IDXMA1)*RATMA1; 
DNO  =  RAT0MA*F509A(IDXMA2)+F509A(IDXMA1)*RATMA1; 
F513A1  =  RAT0MA*F513A(INDXC2)+F513A(INDXC1)*RATMA1; 
F513A2  =  RAT0MA*F513A(INDXC4)+F513A(INDXC3)*RATMA1; 
CMQ  =  RAT0AA>kF513A2+F513A1*RATAA1; 

F514A1  =  RAT0MA*F514A(INDXC2)+F514A(INDXC1)*RATMA1; 
F514A2  =  RAT0MA*F514A(INDXC4)+F514A(INDXC3)*RATMA1; 


CMAD  =  RAT0AA*F514A2+F514A1*RATAA1; 

F601A1  =  RATOMA*F601A(INDXD2)+F601A(INDXD1)*RATMA1; 
F601A2  =  RATQMA*F601A(INDXD4)+F601A(INDXD3)*RATMA1; 
F601A3  =  RATQMA*F601A(INDXD6)+F601A(INDXD5)*RATMA1; 
F601A4  =  RATQMA*F601A(INDXD8)+F601A(INDXD7)*RATMA1; 
F601B1  =  RATQAA*F601A2+F601A1*RATAA1; 

F601B2  =  RATOAA*F601A4+F601A3*RATAA1; 

F601  =  RAT0BA*F601B2+F601B1*RATBA1; 

CNBl  =  F601*sign(BTA); 

F602A1  =  RAT0MC*F602A(INDXE2)+F602A(INDXE1)*RATMC1; 
F602A2  =  RATQMC*F602A(INDXE4)+F602A(INDXE3)*RATMC1; 
F602A3  =  RATQMC*F602A(INDXE6)+F602A(INDXE5)*RATMC1; 
F602A4  =  RATQMC*F602A(INDXE8)+F602A(INDXE7)*RATMC1; 
F602B1  =  RATQAB*F602A2+F602A1*RATAB1; 

F602B2  =  RAT0AB*F602A4+F602A3*RATAB1; 

CNB2  =  RATOBB*F602B2+F602B1*RATBB1; 

F611A1  =  RAT0MA*F611A(INDXC2)+F611A(INDXC1)*RATMA1; 
F611A2  =  RAT0MA*F611A(INDXC4)+F611A(INDXC3)*RATMA1; 
CNP  =  RAT0AA*F611A2+F611A1*RATAA1; 

F612A1  =  RATQMA*F612A(INDXC2)+F612A(INDXC1)*RATMA1; 
F612A2  =  RATQMA*F612A(INDXC4)+F612A(INDXC3)*RATMA1; 
CNR  =  RAT0AA*F612A2+F612A1*RATAA1; 

F613A1  =  RAT0MA*F613A(INDXC2)+F613A(INDXC1)*RATMA1; 
F613A2  =  RAT0MA*F613A(INDXC4)+F613A(INDXC3)*RATMA1; 
CNDA  =  RAT0AA*F613A2+F613A1*RATAA1; 

F614A1  =  RAT0MA*F614A(INDXC2)+F614A(INDXC1)*RATMA1; 
F614A2  =  RAT0MA*F614A(INDXC4)+F614A(INDXC3)*RATMA1; 
CNDD  =  RATQAA*F614A2+F614A1*RATAA1; 

F615A1  =  RAT0MA*F615A(INDXH2)+F615A(INDXH1)*RATMA1; 
F615A2  =  RATDMA*F615A(INDXH4)+F615A(INDXH3)*RATMA1; 
F615A3  =  RAT0MA*F615A(INDXH6)+F615A(INDXH5)*RATMA1; 
F615A4  =  RAT0MA*F615A(INDXH8)+F615A(INDXH7)*RATMA1; 
F615B1  =  RAT0MA*F615A2+F615A1*RATMA1; 
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F615B2  =  RAT0MA*F615A4+F615A3*RATMA1; 

CNDR  =  RATDBD*F615B2+F615B1*RATBD1; 

F617A1  =  RAT0MA*F617A(INDXC2)+F617A(INDXC1)*RATMA1; 
F617A2  =  RAT0MA*F617A(INDXC4)+F617A(INDXC3)*RATMA1; 
DRNK  =  RAT0AA*F617A2+F617A1*RATAA1; 

F618A1  =  RAT0MA*F618A(INDXG2)+F618A(INDXG1)*RATMA1; 
F618A2  =  RATQMA*F618A(INDXG4)+F618A(INDXG3)*RATMA1; 
F618A3  =  RATQMA*F618A(INDXG6)+F618A(INDXG5)*RATMA1; 
F618A4  =  RAT0MA*F618A(INDXG8)+F618A(INDXG7)*RATMA1; 
F618B1  =  RAT0AA*F618A2+F618A1*RATAA1; 

F618B2  =  RAT0AA*F618A4+F618A3*RATAA1; 

DCNSB  =  RAT0BC*F618B2+F618B1*RATBC1; 

IVh  LOOK  UP  CD  AS  A  FUNCTION  OF  CLFT  BASIC  (FlOl)  III 
XDXCLl  =  (CLFTB+1.0)/0.2+1.01; 

IDXCLl  =  floor(XDXCLl): 
if  IDXCLl  <  1,  IDXCLl  =1;  end 
if  IDXCLl  >  13,  IDXCLl  =  13;  end 
INDXBl  =  IDXMA1+10*(IDXCL1-1); 

INDXB2  =  INDXBl+1; 

INDXB3  =  INDXBl+10; 

INDXB4  =  INDXB2+10; 

RATOCL  =  (CLFTB-CLX(IDXCLl))/0.2; 

RATCLl  =  l.O-RATOCL; 

F201A1  =  RAT0MA*F201A(INDXB2)+F201A(INDXB1)*RATMA1; 
F201A2  =  RAT0MA*F201A(INDXB4)+F201A(INDXB3)*RATMA1; 
CDB  =  (RAT0CL*F201A2+F201A1*RATCL1); 

III  CREATE  ’FMCOEF’  VECTOR  III 

FMCOEF(l)  =  CLFTB;  FMC0EF(2)  =  DCLNZ;  FMC0EF(3)  = 

FMCOEF(4)  =  DCLAl;  FMC0EF(5)  =  DCLA2;  FMC0EF(6)  = 

FMC0EF(7)  =  DCDNOZ;  FMC0EF(8)  =  DCDSB;  FMC0EF(9)  = 

FMCOEF(IO)  =  CYB2;  FMCOEF(ll)  =  CYDA;  FMC0EF(12)  = 


DCLO; 

CDB; 

CYBl; 

CYDD; 
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FMC0EF(13)  =  DCYDR; 
FMC0EF(16)  =  CLB2; 
FMC0EF(19)  =  CLDA; 
FMCDEF(22)  =  DRLK; 
FMC0EF(25)  =  DCMNZ; 
FMC0EF(28)  =  DNO; 
FMC0EF(31)  =  CNBl; 
FMCQEF(34)  =  CNR; 
FMCDEF(37)  =  CNDR; 


FMC0EF(14)  =  DRYK; 
FMCDEF(17)  =  CLP; 
FMC0EF(20)  =  CLOD; 
FMC0EF(23)  =  DCLSB; 
FMC0EF(26)  =  DCMO; 
FMC0EF(29)  =  CMQ; 
FMCQEF(32)  =  CNB2; 
FMCQEF(35)  =  CNDA; 
FMCQEF(38)  =  DRNK; 


FMC0EF(15)  =  CLBl; 
FMC0EF(18)  =  CLR; 
FMC0EF(21)  =  DCLDR; 
FMCOEF(24)  =  CMB; 
FMC0EF(27)  =  DNOSB; 
FMCQEF(30)  =  CHAD; 
FMC0EF(33)  =  CNP; 
FMC0EF(36)  =  CNDD; 
FMCDEF(39)  =  DCNSB; 


E.3  F-15  Nonlinear  Propulsion  Model  Listing 

The  following  is  a  listing  of  the  modules  which  comprise  the  nonlinear  propulsion 
model  for  the  F-15. 


•  MATLAB  function:  f25eng 

function  [FPX,FPY,FPZ,DCL,DCM,DCN,TAUL,TAUR,PLAL,PLAR,C0UT1C. . . . 
CQUT2C, FIRST]  =  f 25eng(A, lA) ; 

W/,  INITIALIZATION  SECTION  7.%% 

STEP  =  A(IOOO); 
engdin(STEP)  ; 


7,7,7,  DYNAMICS  (TRIM,  HOLD,  RUN,  AND  LINEARIZATION)  7,7,7, 

[FPX,FPY,FPZ,DCL,DCM,DCN,TAUL,TAUR,PLAL,PLAR,C0UT1C, . . . 
C0UT2C, FIRST]  =  engine(A,IA) ; 

•  MATLAB  function:  engdin 


function  engdin(STEP) ; 


7,7,7,  ROUTINE  TO  READ  IN  ENGINE  DATA  FOR  ENGINE  MODEL  7,7,7, 


124 


III  WRITTEN  AUG.  1,  1978  III 
III  LEE  DUKE  NASA  DFRC  III 

III  COMMON  BLOCKS 
III  ENGVAR 

global  C03RTN  C03RTP  C05F1K  C06LMN  C06LMP  B02LMN  B02LMP 
global  B03RTN  B03RTP  B04F1K  B06F1E  B06F1K 
III  ENGTAB 

global  C03DAT  C03PA  C03NA  C07DAT  B03DAT  B03PA  B03NA  B07DAT. .. 
FIDLAX  FMILAX  FMAXAX 

global  TIDLA  TMILA  TMAXA  FIDLA  FMILA  FMAXA 

HI  =  STEP; 

G060  =  115826.4; 
arg  =  0; 

CC03DAT , C07DAT , B03DAT , B07DAT , TIDLA , FIDLAX , TMILA , FMILAX , TMAXA , . 
FMAXAX]  =  engdat(arg) ; 

FIDLAX  =  FIDLAX/G060; 

FMILAX  =  FMILAX/G060; 

FMAXAX  =  FMAXAX/G060; 

III  DEFINE  CONSTANTS  III 
C05F1C  =  1.0; 

C06LMN  =  0.0; 

C06LMP  =  1.0; 

B02LMN  =  -0.02; 

B02LMP  =  0.02; 

B04F1C  =  1.0; 

B06F1C  =  20.0; 

III  ENTRY  ENGDNl  III 

III  CALCULATIONS  DEPENDENT  ON  FRAME  TIME  III 
C05F1K  =  (C05FlC*HI)/2.0; 


B04F1K  =  (B04FlC*Hl)/2.0; 

B06F1E  =  exp(-B06FlC*Hl); 

B06F1K  =  (l.0-B06FlE)/2.0; 

FIDLA  =  FIDLAX; 

FMILA  =  FMILAX; 

FMAXA  =  FMAXAX; 

•/•/;/.  ADDED,  DUTTON  1/20/94 
C03PA  =  C03DAT(1:4); 

C03NA  =  C03DAT(5:8): 

B03PA  =  B03DAT(1:4); 

B03NA  =  B03DAT(5;8); 

•  MATLAB  function:  engine 

function  [FPX,FPY,FPZ,DCL.DCM,DCN,TAUL,TAUR,PLAL,PLAR, . 
C0UT1C,C0UT2C, FIRST]  *  engine(A,IA) ; 

tVh  ROUTINE  TO  COMPUTE  EFFECTS  OF  ENGINE (S) 

'/,'/'/  WRITTEN  FEBRUARY  1,1983 
Vhl  LEE  DUKE  NASA/DFRF 

Vht  MODIFIED  FOR  DIFFERENT  ANGLE  DEFINITION  JUNE,  1985 

til  JOE  PAHLE  NASA/DFRF 

III  MODIFIED  9/91 

III  JIM  BUFFINGTON  WL/FIGC 

III  VARIABLES  DEFINITIONS: 

III  THRUST ( I, J)  THE  TOTAL  THRUST  OF  THE  I-TH 
III  ENGINE  (LBS) 

III  TLOCAT(I,J)  THE  DISPLACEMENT  OF  THE  I-TH  ENGINE 

III  FROM  THE  AIRCRAFT  REFERENCE 

III  C.G.  (FEET),  WHERE 

III  J=1  CORRESPONDS  TO  THE 

III  THE  X-BODY  AXIS 

III  COORDINATE  OF  THE 


I-TH  ENGINE 

J=2  CORRESPONDS  TO  THE 
THE  Y-BODY  AXIS 
COORDINATE  OF  THE 
I-TH  ENGINE 
J=3  CORRESPONDS  TO  THE 
THE  Z-BODY  AXIS 
COORDINATE  OF  THE 
I-TH  ENGINE 

XYANGL(I)  THE  ANGLE  IN  THE  X-Y  BODY  AXIS 
PLANE  FROM  THE  X-BODY  AXIS  TO 

THE  PROJECTION  OF -THE  I-TH  ENGINE 
AXIS  ONTO  THE  X-Y  PLANE.  (DEG) 

XZANGL(I)  THE  ANGLE  BETWEEN  THE  X-Y  BODY  AXIS 

PLANE  AND  THE  I-TH  ENGINE  AXIS 
PROJECTED  ON  THE  X-Z'  BODY/ENGINE  AXIS 
EIX(I)  THE  MOMENT  OF  INERTIA  ABOUT  THE 

X-ENGINE  AXIS  OF  THE  I-TH  ENGINE 
(SLUG-FT**2) . 

AMSEMG(I)  THE  MASS  OF  THE  ROTATING  MACHINERY 

IN  THE  I-TH  ENGINE  (SLUGS) 

ENGOMG(I)  THE  ROTATIONAL  VELOCITY  OF  THE 

I-TH  ENGINE.  POSITIVE  ROTATION  IS 
MEASURED  USING  THE  R.H.R.  ABOUT  THE 
X-ENGINE  AXIS.  (RAD/SEC) 

TVANXY(I)  THE  ANGLE  IN  THE  X-Y  ENGINE  AXIS 
PLANE  FROM  THE  X-ENGINE  AXIS  TO 

THE  PROJECTION  OF  THE  I-TH  THRUST  VECTOR 
ONTO  THE  X-Y  ENGINE  PLANE.  (DEG) 
TVANXZ(I)  THE  ANGLE  BETWEEN  THE  X-Y  ENGINE  AXIS 

PLANE  AND  THE  I-TH  THRUST  VECTOR 
PROJECTED  ON  THE  X-Z'  ENGINE/THRUST  AXIS 


DXTHRS(I) 


THE  DISTANCE  BETWEEN  THE  C.G.  OF  THE 
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III 

ENGINE  AND  THE  THRUST  POINT  MEASURED 

III 

POSITIVE  IN  THE  NEGATIVE  ENGINE  X-AXIS 

III 

DIRECTION 

III  FPX 

TOTAL 

THRUST  IN  THE  X-BODY  AXIS 

III  FPY 

TOTAL 

THRUST  IN  THE  Y-BODY  AXIS 

III  FPZ 

TOTAL 

THRUST  IN  THE  Z-BODY  AXIS 

III 

DCM 

TOTAL  PITCHING  MOMENT  INCREMENT 

III 

DUE  TO  ENGINES 

III 

DCL 

TOTAL  ROLLING  MOMENT  INCREMENT 

in 

DUE  TO  ENGINES 

III 

DCN 

TOTAL  YAWING  MOMENT  INCREMENT 

III 

DUE  TO  ENGINES 

III  COMMON  BLOCK 

global  THRUST  TLOCAT  XYANGL  XZANGL  TVANXY  TVANXZ  DXTHRS  EIX 
global  AMSENG  ENGOMG 

VH  ASSIGN  A- ARRAY  VAR  NAMES  III 
DGR  =  A(981); 
q  =  A(862); 

P  =  A(861); 

R  =  A(863); 

III  CALL  USER  ENGINE  MODEL  INTERFACE  ROUTINE  III 
I  FIRST  was  originally  set  to  Logical  1  stored  in  data  IA(502) 
CTAUL,TAUR,PLAL,PLAR,C0UT1C,C0UT2C, FIRST]  =  uengin(A,IA) ; 

A(1419)  =  TAUL;  A(1420)  =  TAUR;  A(l43l)  =  PLAL; 

A(667)  =  COUTIC;  A(668)  =  C0UT2C:  IA(502)=FIRST; 

A (1432)  =  PLAR; 

III  COMPUTE  COMPONENTS  OF  THRUST  III 
for  I  =  1:4, 
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ANGLXZ  =  XZAMGL(I)/DGR; 

ANGLXY  =  XYANGL(I)/DGR; 

CSANXZ(I)  =  cos (ANGLXZ); 

SMANXZ(I)  =  sin(ANGLXZ); 

CSANXY(I)  =  cos(AMGLXY); 

SNANXY(I)  =  sin(AMGLXY); 

TANGXY  =  TVAMXY(I)/DGR; 

TAMGXZ  =  TVANXZ(I)/DGR; 

CSTAXY  =  cos (TANGXY); 

SNTAXY  =  sin (TANGXY); 

CSTAXZ  =  cos(TANGXZ); 

SNTAXZ  =  sin(TANGXZ); 

XTHRSI  =  THRUST(I)*CSTAXZ*CSTAXY; 

YTHRSI  =  THRUST(I)*CSTAXZ*SNTAXY; 

ZTHRSI  =  -THRUST(I)*SNTAXZ; 

XTHRST(I)  =  XTHRSI*CSANXZ (I) *CSANXY(I) -YTHRSI*. . . 

SNANXY(I)+ZTHRSI*SNANXZ(I)*CSANXY(I) ; 
YTHRST(I)  =  XTHRSI*CSANXZ(I)*SNANXY(I)+YTHRSI*. . . 

CS ANXY (I ) +ZTHRSI*SNANXZ ( I ) *CS ANXY (I ) ; 
ZTHRST(I)  =  -XTHRSI*SNANXZ(I)+ZTHRSI*CSANXZ(I); 

end 


III  COMPUTE  TOTAL  X-AXIS  AND  Z-AXIS  THRUST  HI 
FPX  =  0.0; 

FPY  =  0.0; 

FPZ  =  0.0; 
for  I  =  1:4, 

FPX  =  FPX+XTHRST(I); 

FPY  =  FPY+YTHRST(I); 

FPZ  =  FPZ+ZTHRST(I) ; 

end 


III  COMPUTE  ROTATIONAL  EFFECTS  (TORQUE)  DUE  TO  ENGINE  OFFSET  HI 


yxi  FROM  CENTERLINE,  AND  THRUST  VECTORING  III 

TORQUE  =  zeros (3,1); 
for  I  =  1:4, 

TVECLCd,!)  =  TL0CAT(I,1)-DXTHRS(I)*CSANXZ(I)*CSANXY(I); 
TVECLC(I,2)  =  TL0CAT(I,2)-DXTHRS(I)*CSANXZ(I)*SNANXY(I); 
TVECLCd, 3)  =  TL0CATd,3)+DXTHRSd)*SNANXZd)  ; 

TORQUEd)  =  T0RQUE(l)+ZTHRSTd)*TVECLC(I,2)-.  .  . 
YTHRST(I)*TVECLCd,3)  ; 

T0RQUE(2)  =  T0RQUE(2)+XTHRSTd)*TVECLC(l,3)-.  .  . 

ZTHRST (1) ♦TVECLC (I , 1) ; 

T0RqUE(3)  =  TORQUE(3)+YTHRSTd)*TVECLCd,l)-.  .  . 

XTHRST ( I ) *TVECLC (1,2); 

end 

III  COMPUTE  ANGULAR  MOMENTUM  OF  ENGINES  III 
for  I  =  1:4, 

OMEGX  =  ENGOMGd)*CSANXZd)*CSANXYd)  ; 

OMEGY  =  ENGOMGd)*CSANXZd)*SNANXYd)  ; 

OMEGZ  =  -ENGOMGd)*SNANXZd); 

EIXE  =  ElXd); 

EX  =  TLOCATd,!); 

EY  =  TL0CAT(I,2); 

EZ  =  TLOCATd, 3); 

III  COMPUTE  THE  INERTIA  TENSOR  OF  THE  I-TH  ENGINE  III 
III  ABOUT  THE  AIRCRAFT  REFERENCE  C.G.  Ill 

AINRTA(1,1)  =  EIXE*CSANXZd)*CSANXZd)*.  .  . 

CSANXYd)*CSANXYd)  ; 

AINRTA(1,2)  =  EIXE*CSANXZd)*CSANXZd)*.  .  . 

SNANXYd)*CSANXYd)  ; 

AINRTA(1,3)  =  -EIXE*CSANXZd)*SNANXZd)*.  .  . 

CSANXYd); 

AINRTA(2,1)  =  AINRTA(1,2); 

AINRTA(2,2)  =  EIXE*SNANXYd)*SNANXYd)* .  .  . 
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AINRTA(2,3) 

AINRTA(3,1) 

AINRTA(3,2) 

AINRTA(3,3) 


CSANXZ(I)*CSANXZ(I); 
-EIXE*SNANXZ(I)*CSANXZ(I)*SNANXY(I) ; 
AIMRTA(1,3); 

AIMRTA(2,3): 

EIXE*SNANXZ(I)*SNANXZ(I) ; 


for  J  =  1:3, 

HEMG IN ( I , J )  =  QMEGX+AINRTA ( J , 1) + DMEGY* . . . 

AINRTA(J,2)+DMEGZ*AIMRTA(J,3) ; 


end 


end 


COMPUTE  GYROSCOPIC  EFFECTS  VLt 
GYRO  =  zeros (3,1); 
for  I  =  1:4, 

GYRO(l)  =  GYR0(1)+Q*HENGIN(I,3)-R*HENGIN(I,2) ; 
GYR0(2)  =  GYR0(2)+R*HENGIN(I,1)-P*HENGIN(I,3) ; 
GYR0(3)  =  GYR0(3)+P:t:HENGIN(I,2)-Q*HENGIN(I,l)  ; 

end 


'/’/,•/  COMPUTE  TOTAL  MOMENT  INCREMENT  DUE  TO  ENGINES  Vhl 

DCL  =  GYRO(l)+TORqUE(l); 

DCM  =  GYR0(2)+T0RQUE(2); 

DCN  =  GYR0(3)+T0RQUE(3); 


•  MATLAB  function:  uengin 


function  [TAUL,TAUR,PLAL,PLAR,C0UT1C,C0UT2C, FIRST]  =  uengin(A,IA) ; 

IVL  ROUTINE  TO  MODEL  THE  ENGINE  VtX 

FIRST  =  IA(502);  VH  Logical  TRUE  III 

III  COMMON  BLOCK 

global  THRUST  TLOCAT  XYANGL  XZANGL  TVANXY  TVANXZ  DXTHRS 
global  EIX  AMSENG  ENGOMG 
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if  FIRST==1 


FIRST  =  0;  III  Logical  FALSE  III 


THRUST  =  zeros (4,1); 
TLOCAT  =  zeros (4, 3); 
XYANGL  =  zeros (4,1); 


XZANGL  =  zeros (4,1); 


TVANXY  =  zeros (4,1); 
TVANXZ  =  zeros (4,1); 
DXTHRS  =  zeros (4,1); 
ENGOMG  =  zeros (4,1); 
EIX  =  zeros(4,l) ; 
III  ENGINES:  2  ENGINES  10  FT. 

TLOCATd,!)  =  -10.0; 


III  ADDED,  DUTTON  1/20/94 
III  ADDED,  DUTTON  1/20/94 
III  ADDED,  DUTTON  1/20/94 
III  ADDED,  DUTTON  1/20/94 
III  ADDED,  DUTTON  1/20/94 
BEHIND  CG,  4  FT  OFF  CENTERLINE  III 


TLOCATd, 2)  =  -4.0; 
TL0CAT(2,1)  =  -10.0; 
TL0CAT(2,2)  =  4.0; 


end 


III  CALL  ROUTINE  TO  MODEL  ENGINE  III 
[TAUL,TAUR,PLAL,PLAR,C0UT1C,C0UT2C]  =  engmdl(A,IA) ; 


•  MATLAB  function:  engmdl 

function  [TAUL,TAUR,PLAL,PLAR,C0UT1C,C0UT2C]  =  engmdl(A,IA) ; 

III  ROUTINE  TO  IMPLEMENT  FIRST  ORDER  ENGINE  MODEL  III 
III  WRITTEN  AUG.  2,  1978  III 
III  LEE  DUKE  NASA/DFRC  III 

III  COMMON  BLOCKS 
III  ENGVAR 

global  C03RTN  C03RTP  C05F1K  C06LMN  C06LMP  B02LMN  B02LMP 
global  B03RTN  B03RTP  B04F1K  B06F1E  B06F1K 
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III  TLUENG 

global  TIDL  TMIL  TMIN  TMAX  FIDL  FMIL  FMIN  FMAX 
III  CONPOS 

global  DAP  FLAT  DATRIM  DEP  FLQN  DETRIM  DRP  FPED  DRTRIM  DSBP 
global  DFP  THSL  THSR 
III  ENGSTF 

global  THRUST  TLDCAT  XYANGL  XZAWGL  TVANXY  TVANXZ  DXTHRS 
global  EIX  AMSEMG  EMGOMG 

III  ASSIGN  A,IA-ARRAY  VAR  NAMES  HI 
IMQDE  =  IA(501);  PLAPL  =  A(1416); 

PLAPR  =  A(1417);  PLASYM  =  A(1418); 


III  ADDED,  DUTTON  4/7/94  III 
CDUTIC  =  A(667);  C0UT2C  =  A(668); 

III  ASSIGN  DATA  VALUES  HI 
CORMIN  =  20.0;  CORDIS  =  63.0; 

AUGMIN  =  83.1;  AUGDIS  =  44.0; 


FMSSIC 

= 

0.0; 

B04IN2 

= 

0.0;  B040T2 

= 

0.0; 

B06IN2 

= 

0.0;  B060T2 

= 

0.0; 

B14IN2 

= 

0.0;  B140T2 

= 

0.0; 

B16IN2 

= 

0.0;  B160T2 

= 

0.0; 

C05IN2 

= 

0.0;  C050T2 

= 

0.0; 

C060T2 

= 

0.0;  C070T2 

= 

0.0; 

C15IN2 

= 

0.0;  C150T2 

= 

0.0; 

C160T2 

= 

0.0;  C170T2 

= 

0.0; 

ENGLU  =  1.02; 
ENGRU  =  .97; 
if  PLASYM  >  le-10 
PLAPL  =  PLASYM; 
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PLAPR  =  PLASYM; 
ENGLU  =1.0; 
ENGRU  =1.0; 
end 


'/;/;/  FIND  THRUST  AND  FUEL  FLOW  FOR  PRESENT  lit 
III  MACH/ ALTITUDE  CONDITION  III 
H  =  A(713); 

AMCH  =  A(825); 
engtlu(H,AMCH) ; 

III  ENGINE  COMMAND  INPUTS  III 
III  RIGHT  ENGINE  CALCULATIONS  III 

III  DETERMINE  INPUTS  TO  CORE  AND  AUGMENTOR  MODELS  III 
CINOl  =  (PLAPR-CORMIN) /CORDIS; 

[CINOl]  =  f limit (CINOl,  0.0,  1.0); 

BINOl  =  (PLAPR-AUGMIN)/AUGDIS; 

[BINOl]  =  f limit (BINOl,  0.0,  1.0); 

III  WAIT  FOR  CORE  TO  SPOOL  UP  TO  99%  BEFORE  BURNER  IS  ACTIVE  %%% 
if  COUTIC  <0.99 
BINOl  =0.0; 

B040T2  =0.0; 

B05TIM  =0.0; 
end 


III  FIRST  ORDER  CORE  ENGINE  MODEL  III 
if  IM0DE==5 

C02JN1  =  CIN01-C060T2; 

C03IN1  =  C02JN1; 

[C030T1]  =  f limit (C03IN1, C03RTN, C03RTP) ; 
C04IN1  =  C030T1; 

C04MLT  =  C070T2; 
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C040T1  =  C04IN1*C04MLT; 

C05IN1  =  C04QT1; 

C050T1  =  C05QT2+C05F1K*(C05IN1+C05IN2) ; 
else 

C050T1  =  CINOl; 
end 

C06IN1  =  C050T1; 

CC060T1]  =  f limit (C06IN1, C06LMN, C06LMP) ; 
C07IN1  =  C05QT1; 

[C07QT1]  =  c07sdl(C07IMl); 
if  abs(C07QT2)  >=  le-10 
TAUR  =  1/C070T2; 
else 

TAUR  =  lelO; 
end 

COUTIC  =  C060T1; 

III  SIMPLE  AUGMENTOR  MODEL  III 
if  IM0DE==5 

BOIJNI  =  BIN01-B050T2; 

B02IN1  =  BOIJNI; 

B020T1  =  B02IN1; 

B03IN1  =  B020T1; 

[B030T1]  =  flimit(B03INl,B03RTN,B03RTP) ; 
B04IN1  =  B030T1; 

B040T1  =  B040T2+B04F1K*(B04IN1+B04IN2) ; 
else 

B040T1  =  BINOl; 
end 

CB040T1]  =  f limit (B040T1, 0.0, 1.0); 

B050T1  =  B040T1; 

B06IN1  =  B050T1; 
if  IM0DE==5 


B060T1  =  B060T2*B06F1E+B06F1K*(B06IN1+B06IN2); 
else 

B06QT1  =  B05QT1; 
end 

B07IM1  =  B060T1; 

[B070T1]  =  b07sdl(B07INl); 

B08IN1  =  B070T1; 

[B08QT1]  =  b08f cn(B08IMl) ; 

BOUTIC  =  B08QT1; 

B0UT2C  =  B070T1; 

VH  OUTPUTS  lit 

CQROTR  =  COUTIC; 

AUGOTR  =  BOUTIC; 

AUGPLR  =  B0UT2C; 

III  UPDATE  FOR  NEXT  FRAME  III 
C050T2  =  C050T1; 

C05IN2  =  C05IN1; 

C060T2  =  C060T1; 

C070T2  =  C070T1; 

B040T2  =  B040T1; 

B04IN2  =  B04IN1; 

B050T2  =  B050T1; 

B060T2  =  B060T1; 

B06IN2  =  B06IN1; 

III  LEFT  ENGINE  CALCULATIONS  III 

III  DETERMINE  INPUTS  TO  CORE  AND  AUGMENTOR  MODELS  III 
CINll  =  (PLAPL-CORMIN) /CORDIS; 

[CINll]  =  f limit (CINll, 0.0, 1.0); 

BINll  =  (PLAPL-AUGMIN)/AUGDIS; 

[BINll]  =  flimit(BINll,0.0,1.0); 
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VH  CORE  TO  AUGMENTOR  SWITCHING  LOGIC  lit 

III  WAIT  FOR  CORE  TO  SPOOL  UP  TO  99%  BEFORE  BURNER  IS  ACTIVE  III 
if  C0UT2C  <0.99 
BINll  =  0.0; 

B140T2  =0.0; 

B15TIM  =0.0; 
end 


III  FIRST  ORDER  CORE  ENGINE  MODEL  III 
if  IM0DE==5 

C12JN1  =  CIN11-C16QT2; 

C13IN1  =  C12JN1; 

CC130T1]  =  flimit(C13INl,C03RTM,C03RTP); 
C14IN1  =  C130T1; 

C14MLT  =  C170T2; 

C140T1  =  C14INl>t<C14MLT; 

C15IN1  =  C140T1; 

C150T1  =  C150T2+C05F1K*(C15IN1+C15IN2) ; 
else 

C150T1  =  CINll; 
end 

C16IN1  =  C150T1; 

[C160T1]  =  f limit (C16IN1, C06LMN, C06LMP) ; 
C17IN1  =  C150T1; 

[C170T1]  =  c07sdl(C17INl); 
if  abs(C170T2)  >=  le-10 
TAUL  =  1/C170T2; 
else 

TAUL  =  lelO; 
end 

C0UT2C  =  C160T1; 
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III  SIMPLE  AUGMENTOR  MODEL  TLl 
if  IM0DE==5 

BllJMl  =  BIN11-B15QT2; 

B12IN1  =  BllJNl; 

B120T1  =  B12IN1; 

B13IM1  =  B120T1; 

[B13QT1]  =  f limit (B13IN1, B03RTN, B03RTP) ; 

B14IN1  =  B130T1; 

B14QT1  =  B140T2+B04F1K*(B14IN1+B14IN2) ; 
else 

B140T1  =  BINll; 
end 

[B140T1]  =  f limit (B140T1, 0. 0, 1. 0) ; 

B15QT1  =  B140T1; 

B16IN1  =  B150T1; 
if  IM0DE==5 

B160T1  =  B160T2*B06F1E+B06F1K*(B16IN1+B16IN2); 
else 

B160T1  =  B150T1; 
end 

B17IN1  =  B160T1; 

[B170T1]  =  b07sdl(B17INl); 

B18IN1  =  B170T1; 

[B180T1]  =  b08fcn(B18INl); 

B0UT3C  =  B180T1; 

B0UT4C  =  B170T1; 

III  OUTPUTS  III 

COROTL  =  C0UT2C; 

AUGOTL  =  B0UT3C; 

AUGPLL  =  B0UT4C; 


III  UPDATE  FOR  NEXT  FRAME  III 
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C150T2  =  C150T1; 

C15IN2  =  C15IN1; 

C160T2  =  C160T1: 

C170T2  =  C170T1; 

B14DT2  =  B14QT1; 

B14IM2  =  B14IN1; 

B150T2  =  B150T1; 

B160T2  =  B160T1; 

B16IN2  =  B16IN1; 

tit  COMPUTE  ENGINE  MODEL  OUTPUTS  ttt 
if  COROTR  <0.99 

PLAR  =  COROTR*CORDIS+CORMIN; 

FLOWR  =  (FMIL-FIDL)*COROTR+FIDL; 

THSR  =  (TMIL-TIDL)*COROTR+TIDL; 
else 

PLAR  =  COROTR*CORDIS+AUGOTR*AUGDIS+CORMIN; 

FLOWR  =  (FMAX-FMIL)*AUGOTR+(FMIL-FIDL)*COROTR+FIDL; 
THSR  =  (TMAX-TMIL)*AUGOTR+(TMIL-TIDL)*COROTR+TIDL; 
end 

if  COROTR  <0.99 

PLAL  =  COROTL*CORDIS+CORMIN; 

FLOWL  =  (FMIL-FIDL)*COROTL+FIDL; 

THSL  =  (TMIL-TIDL)*COROTL+TIDL; 
else 

PLAL  =  COROTL*CORDIS+AUGOTL*AUGDIS+CORMIN; 

FLOWL  =  (FMAX-FMIL)*AUGOTL+(FMIL-FIDL)*COROTL+FIDL; 
THSL  =  (TMAX-TMID^AUGOTL+CTMIL-TIDD+COROTL+TIDL; 
end 


ttt  STORE  IN  COMMON  ttt 

THRUST(l)  =  ENGRU*THSR; 
THRUST (2)  =  ENGLU*THSL; 


139 


E.4  Atmospheric  Model  Listing 

The  following  is  a  listing  of  the  modules  which  comprise  the  atmospheric  model  used 
in  the  nonlinear  F-15  simulation. 


•  MATLAB  function:  atmos 


function  [AMCH,RHQ ,QBAR,G]  =  atmos (A); 

lit  ROUTINE  TO  CALCULATE  AIR  DATA  PARAMETERS  HI 
III  WRITTEN  MAY  13,  1981  III 
III  LEE  DUKE  NASA/DFRC  III 

III  ASSIGN  A- ARRAY  VAR  NAMES  HI 
V  =  A<829); 

H  =  A(713); 

III  INVOKE  ATMOSPHERIC  MODEL  IH 

[An,RHO,G,PA,TMPR,VMU]  =  altfn(H) ; 

HI  COMPUTE  MACH  HI 
AMCH  =  V/An; 

III  COMPUTE  DYNAMIC  PRESSURE  HI 
QBAR  =  RH0*V''2/2.0; 


•  MATLAB  function:  altfn 


function  [A,RHO,G,PA,TMPR,VMU]  =  altfn (H) ; 


III  ROUTINE  PROVIDES  BY  TABLE  LOOK  UP  HI 
III  VELOCITY  OF  SOUND  A  -  (FT/SEC)  HI 
III  ACCELERATION  DUE  TO  GRAVITY  G  -  (FT/SEC>t=*2)  III 
III  AIR  DENSITY  RHO  -  (SLUGS/FT**3)  HI 
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11%  AMBIENT  STATIC  PRESSURE  PA  -  (PSF)  111 

11%  AMBIENT  AIR  TEMPERATURE  TMPR  -  (DEG  RANKIN)  111 

%X%  AS  A  FUNCTION  OF  ALTITUDE,  H  (FT) 

%'/,%  BREAK  POINTS  ON  H  ARE  3280.84  FT  (1000. 0  M) 

%XX<,  The  data  array  is  of  size  81. 

*/////»  The  first  data  XX(l)  corresponds  to  SEA  LEVEL  to  999m 
'/o’/o'/o  The  second  data  XX(2)  corresponds  to  1000m  to  1999m  etc... 
y////o  Equal  to  or  above  80,000m  the  last  array  element  is  XX(81) . 


'%"%%  WRITTEN  11/9/72  A  MYERS  NASA  FRC 

'/'///  TMPR  ADDED  1  MAY  78  L  SCHILLING  NASA/DFRC 

AA  =  [1116.45,1103.79,1090.98,1078.03,1064.92,1051.66,1038.23. 
1024.63,1010.84,  996.88,  982.72,  968.07,  968.07,  968.07. 
968.07,  968.07,  968.07,  968.07,  968.07,  968.07,  968.07.. 

970.15,  972.37,  974.57,  976.77,  978.97,  981.15,  983.34.. 

985.52,  987.69,  989.86,  992.02,  994.28,  999.56,1005.54.. 

1011.48,1017.39,1023.25,1029.08,1034.88,1040.65,1046.38. 
1052 .07,1057.74, 1063 . 37 , 1068 . 97 , 1074 . 54 , 1082 . 02 , 1082 . 02 . 
1082 . 02 , 1082 . 02 , 1082 . 02 , 1082 . 02 , 1079 . 77 , 1075 .82,1071.86. 
1067 . 89 , 1063 . 90 , 1059 . 90 , 1055 .89,1051.86, 1047 . 81 , 1042 . 09 . 
1033.92,1025.68,1017.39,1009.02,1000.59,  992.08,  983.51. 
974.87,  966.15,  957.35,  948.47,  939.52,  930.48,  921.36.. 

912.14,  902.82,  893.44,  883.99]; 

RHOA  =  [.23769D-2, .21571D-2,.19531D-2,.17642D-2, .15898D-2. . . 

. 14289D-2 , . 12808D-2 , . 11448D-2 , . 10202D-2 , . 90625D-3 , . 80234D-3 . 
.70783D-3, .60526D-3, .51729D-3,.44212D-3, .37788D-3, .32301D-3. 
. 27611D-3 , . 23604D-3 , . 20179D-3 , . 1725 lD-3 , . 14691D-3 , . 12517D-3 . 
.10673D-3, .91075D-4, .77776D-4,.66470D-4, .56848D-4, .48655D-4. 
.41674D-4, .35721D-4, .30642D-4, .26301D-4, .22455D-4, . 19185D-4. 
.16422D-4, .14083D-4, .12099D-4,.10413D-4, .89773D-5, .77529D-5 . 
. 67065D-5 , . 58 109D-5 , . 50427D-5 , . 43830D-5 , . 38153D-5 , . 33259D-5 . 


. 29037D-5 , . 25548D-5 , . 22562D-5 , . 19925D-5 , . 17597D-5 , . 15541D-5 . 
. 13782D-5 , . 12251D-5 , . 10880D-5 , . 96554D-6 , . 85609D-6 , . 75839D-6 . 
. 67123D-6 , . 59358D-6 , . 52443D-6 , . 46434D-6 , . 41236D-6 , . 36554D-6 . 
. 32335D-6 , . 28548D-6 , . 25154D-6 , . 22118D-6 , . 19403D-6 , . 16985D-6 . 
. 14823D-6 , . 12921D-6 , . 11228D-6 , . 97309D-7 , . 841 13D-7 , . 72490D-7 . 
.  62284D-7 , . 53359D-7 , . 45578D-7 , . 38787D-7] ; 

GA  =  [32.174,32.164,32.154,32.144,32.134,32.123,32.114,32.103.. 
32.093,32.083,32.073,32.063,32.053,32.043,32.033,32.023. . 
32.013,32.003,31.992,31.983,31.972,31.963,31.952,31.943. . 

31.932.31.923.31.912.31.903.31.892.31.883.31.872.31.863. . 

31.852.31.843.31.833.31.823.31.813.31.803.31.793.31.783. . 

31.773.31.763.31.753.31.743.31.733.31.723.31.713.31.703. . 

31.694.31.684.31.674.31.664.31.654.31.644.31.634.31.624. . 

31.615.31. 605.31 . 595.31 . 585 .31.575.31.565.31. 555 .31.546.. 

31.536.31.526.31.516.31.506.31.496.31.487.31.477.31.467. . 

31.457.31.448.31.438.31.427.31.417.31.407.31.398. .  . . 
31.388,31.378] ; 

PAA  =  [.211622D+4,.187711D+4,.166042D+4,.146451D+4,.128781D+4.. 

. 112882D+4, .986161D+3, .858503D+3, .744600D+3, .643286D+3. . 
. 55346 lD+3, .474097D+3,. 405 167D+3, .3462720+3, .2959530+3. . 
. 2529610+3 , . 2 162230+3 , . 1848300+3 , . 1580030+3 , . 1350760+3 . . 
. 1154820+3 , . 9876590+2 , . 8453380+2 , . 7240700+2 , . 6206620+2 . . 
. 5324160+2 , .4570510+2 , . 3926400+2 , . 3375510+2 , . 2903960+2 . . 
. 2500050+2 , .2153830+2 , . 1856850+2 , . 1602560+2 , . 1385570+2 . . 
. 1200060+2 , . 1041180+2 , . 9048570+1 , . 7876710+1 , . 6867580+1 . . 
. 5997120+1 , . 5245020+1 , .4594120+1 , . 4029870+1 , . 3540020+1 . . 
. 3114040+1 , . 2743110+1 , . 2419600+1 , . 2136490+1 , . 1886720+1 . . 
. 1666220+1 , . 1471550+1 , . 1299670+1 , . 1147740+1 , . 1012770+1 . . 
. 8928890+0 , . 7864880+0 , . 6921380+0 , . 6085470+0 , . 5345460+0 . . 
. 4691020+0 , .4112710+0 , . 3601830+0 , . 3148570+0 , . 2746530+0 . . 
. 2390620+0 , . 2076160+0 , . 1798880+0 , . 1554890+0 , . 1340660+0 . . 
. 1152980+0 , . 9889320-1 , . 8458890-1 , . 7214700-1 , . 6135260-1 . . 
.5201170-1,. 4395330-1,. 3701990-1,. 3107220-1,. . . 
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.259855D-1, .216500D-1] ; 

TMPRA  =  [518.670,  506.972,  495.277,  483.586,  471.899... 

460.217,  448.537,  436.860,  425.187,  413.159... 
401.854,  390.193,  389.970,  389.970,  389.970... 
389.970,  389.970,  389.970,  389.970,  389.970... 
389.970,  391.646,  393.433,  395.221,  397.008... 
398.794,  400.579,  402.365,  404.149,  405.932... 
407.716,  409.500,  411.282,  415.751,  420.737... 
425.723,  430.708,  435.690,  440.672,  445.651... 
450.630,  455.605,  460.580,  465.554,  470.525... 
475.495,  480.465,  485.431,  487.170,  487.170... 
487.170,  487.170,  487.170,  485.149,  481.608... 
478.069,  474.530,  470.993,  467.458,  463.923... 
460.390,  456.858,  451.883,  444.821,  437.764... 
430.708,  423.653,  416.603,  409.552,  402.505... 
395.460,  388.418,  381.377,  374.339,  367.303... 
360.270,  353.232,  346.212,  339.174,  332.154,  325.170] 

lit  DEFINE  THE  BREAK  POINT 
FTBP  =  3280.84; 
lit  INTERPOLATION  OF  THE  DATA 
if  H  <=  0.0 

A  =  AA(1); 

RHO  =  RHOA(l); 

G  =  GA(1); 

PA  =  PAA(l); 

TMPR  =  TMPRA(l) ; 

III  EQUATION  FOR  VISCOSITY  (VMU)  ACCURATE  TO  90  KM  ONLY  III 
VMU  =  ((7.3025D-07)*(TMPR)~(1.5))/(TMPR+198.72) : 
elseif  H  >=  262467.2 
A  =  AA(81); 

RHO  =RH0A(81); 

G  =  GA(81); 
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PA  =  PAA(81); 

TMPR  =  TMPRA(81) ; 

VISCOSITY  above  90  KM  SET  AT  90  KM  VALUE  lit 
VMU  =  8.32627D-06; 
else 

RBP  =  H/FTBP ; 

11  =  floor(RBP)+l; 

12  =  Il+l; 

R  =  RBP-Il+l; 

A  =  R*(AA(I2)-AA(I1))+AA(I1); 

RHO  =  R*(RHQA(I2)-RH0A(I1))+RH0A(I1) ; 

G  =  R*(GA(I2)-GA(I1))+GA(I1); 

PA  =  R*(PAA(I2)-PAA(I1))+PAA(I1); 

TMPR  =  R*(TMPRA(I2)-TMPRA(I1))+TMPRA(I1) ; 

III  EQUATION  FOR  VISCOSITY  (VMU)  IS  ACCURATE  UP  TO  90  KM  ONLY 
VMU  =  ((7.3025D-07)*(TMPR)~(1.5))/(TMPR+198.72) ; 


end 


Appendix  F 

SAMPLE  SIMULATION  COMMAND  LISTING 

F.l  Nonlinear  Open-Loop  Simulation 

The  following  summarizes  the  flow  of  commands  to  execute  the  nonlinear  open-loop 
F-15  simulation: 

•  Start  MATLAB. 

•  Using  an  editor,  enter  the  desired  trim  values  for  Xo  and  Uo  in  trimmod.m. 

•  Type  f  251oad  in  MATLAB. 

•  Type  f  25sim  in  MATLAB  to  load  the  SIMULINK  environment. 

•  Select  Simulation  menu,  select  Parameters,  and  set  desired  simulation  parame¬ 
ters.  Note:  The  open-loop  simulations  in  this  report  used  the  Euler  method  of 
numerical  integration. 

•  Select  Simulation  menu  and  select  Start. 

•  If  you  desire  to  monitor  simulation’s  progress,  double-click  on  the  Clock  icon 
and  a  running  count  of  the  time  will  be  shown. 

•  Once  the  simulation  is  complete,  the  selected  parameters  are  stored  in  the 
MATLAB  workspace. 

F.2  Nonlinear  Closed-Loop  Simulation 

The  following  summarizes  the  flow  of  commands  to  execute  the  nonlinear  closed-loop 
F-15  simulation: 


•  Start  MATLAB. 
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•  Using  an  editor,  enter  the  desired  trim  values  for  Xq  and  Uq  in  trimmod.m. 

•  Type  f  251oad  in  MATLAB. 

•  Enter  the  values  for  the  TECS  controller  gains.  If  either  FPl  or  FP2  is  be¬ 
ing  used,  these  gains  can  be  loaded  by  typing  load  f25sd539_cx  or  load 
f25sd497_cx  respectively. 

•  Type  f  25simtecs  in  MATLAB  to  load  the  SIMULINK  environment. 

•  Select  Simulation  menu,  select  Parameters,  and  set  desired  simulation  param¬ 
eters.  Note:  The  closed-loop  simulations  in  this  report  used  the  Euler  method 
of  numerical  integration  and  minimum  time  steps  of  no  less  than  0.01. 

•  Select  Simulation  menu  and  select  Start. 

•  If  you  desire  to  monitor  simulation’s  progress,  double-click  on  the  Clock  icon 
and  a  running  count  of  the  time  will  be  shown. 

•  Once  the  simulation  is  complete,  the  selected  parameters  are  stored  in  the 
MATLAB  workspace. 

•  To  quickly  view  the  output,  type  plotcl  in  MATLAB. 

F.3  Linear  Closed-Loop  Simulation 

The  following  summarizes  the  flow  of  commands  to  execute  the  linearized  closed-loop 

F-15  simulation: 

•  Start  MATLAB. 

•  Load  the  state  space  model.  If  either  FPl  or  FP2  is  being  used,  they  can  be 
loaded  by  typing  trim497ss  or  trim539ss  respectively. 

•  Enter  the  values  for  the  TECS  controller  gains.  If  either  FPl  or  FP2  is  be¬ 
ing  used,  these  gains  can  be  loaded  by  typing  load  f25sd539_cx  or  load 
f  25sd497_cx  respectively. 


Type  f  25simtecslin  in  MATLAB  to  load  the  SIMULINK  environment. 


Select  Simulation  menu,  select  Parameters,  and  set  desired  simulation  param¬ 
eters.  Note:  The  closed-loop  simulations  in  this  report  used  the  Euler  method 
of  numerical  integration  and  minimum  time  steps  of  no  less  than  0.01. 

Select  Simulation  menu  and  select  Start. 

If  you  desire  to  monitor  simulation’s  progress,  double-click  on  the  Clock  icon 
and  a  running  count  of  the  time  will  be  shown. 

Once  the  simulation  is  complete,  the  selected  parameters  are  stored  in  the 
MATLAB  workspace. 


To  quickly  view  the  output,  type  plotlin  in  MATLAB. 


